Thank you @ptrblck ! your advices were very helpful.
I report here the explicit code to make it available for users who will need it.
If anyone notice something that seems wrong please don’t esitate to report it.
class VGG_Custom_Dropout(nn.Module, NetVariables, OrthoInit):
    def __init__(self, params):
      
        self.params = params.copy()
        
        nn.Module.__init__(self)
        NetVariables.__init__(self, self.params)
        OrthoInit.__init__(self)
        
        
        self.ModuleDict = self._make_layers(cfg['VGG16']) #CREATE A DICT OF MODULES WITH IDENTITY ON DROPOUT PLACE
        self.classifier = nn.Linear(512, self.num_classes)
        self.mask_dict = {}


        self.weights_init() #call the orthogonal initial condition
        


    def forward(self, x, Mask_Flag, TrainMode_Flag):
        outs = {} 
        #WE ITERATE OVER THE MODULE DICT SUBSTITUTING IDENTITY WITH THE DROPOUT MASK WHEN THE TrainingMode_Flag TELL US WE ARE IN TRAINING MODE
        
        for key in self.ModuleDict:
            if key.startswith('!'):
                
                if not TrainMode_Flag: #if we are in eval mode the dropout is substitute by a identity layer
                    x = self.ModuleDict[key](x)
                else:
                
                
                    if Mask_Flag==1: #this flag trigger the update of masks 
                        self.mask_dict[key] = torch.distributions.Bernoulli(probs=(1-self.params['dropout_p'])).sample(x.size())

                    x = x * self.mask_dict[key] * 1/(1-self.params['dropout_p']) #dropout layer

                    
            else: #for modules different from dropout regular forward
                x = self.ModuleDict[key](x)
        
        L2 = x
        outs['l2'] = L2
        Out = L2.view(L2.size(0), -1)
        Out = self.classifier(Out)
        outs['out'] = Out
        return outs

    def _make_layers(self, cfg):
        layers = []
        in_channels = 3
        ModuleDict = nn.ModuleDict()
        NumberKey=0
        
        for x in cfg:
            if x == 'M':
                layers += [nn.MaxPool2d(kernel_size=2, stride=2)]
                
            elif x=='A':
                layers += [nn.AvgPool2d(kernel_size=2, stride=2)]
                 
                
            elif x == 'Dp':
                
                ModuleDict[str(NumberKey)] = nn.Sequential(*copy.deepcopy(layers)) 
                NumberKey+=1
                ModuleDict['!'+str(NumberKey)] = nn.Identity()
                NumberKey+=1
                layers = []
            else:
                layers += [nn.Conv2d(in_channels, x, kernel_size=3, padding=1),
                           #nn.ReLU(inplace=True),
                           #nn.LeakyReLU(negative_slope=0.1, inplace=False),
                           nn.Tanh()
                           #,nn.BatchNorm2d(x)   #For now Batch Norm is excluded because it is incompatible with PCNGD, GD, PCNSGD where I forward sample by sample
                           ,nn.GroupNorm(int(x/self.params['group_factor']), x)
                           #,nn.GroupNorm(int(1), x)
                           ]
                in_channels = x
        NumberKey+=1
        ModuleDict[str(NumberKey)] = nn.AvgPool2d(kernel_size=1, stride=1)
        return ModuleDict

In this way I have control on where the mask has to change and when it has to be substitued by an identity layer (when we are in the evaluation mode). Do I miss something? For example should I detach the mask multiplication?
I think it should not make difference because also if the gradient record that operation I’m turning off the nodes, so no signal leave from them. Is it right?