8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors
534 write (iout,*) 'FTORS factor =',ftors
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
556 if ( secstruc(i) .eq. 'H') then
557 ! Helix restraints for this residue
560 phi0(ii) = 45.0D0*deg2rad
561 drange(ii)= 5.0D0*deg2rad
562 phibound(1,i) = phi0(ii)-drange(ii)
563 phibound(2,i) = phi0(ii)+drange(ii)
564 else if (secstruc(i) .eq. 'E') then
565 ! strand restraints for this residue
568 phi0(ii) = 180.0D0*deg2rad
569 drange(ii)= 5.0D0*deg2rad
570 phibound(1,i) = phi0(ii)-drange(ii)
571 phibound(2,i) = phi0(ii)+drange(ii)
573 ! no restraints for this residue
574 ndih_nconstr=ndih_nconstr+1
575 idih_nconstr(ndih_nconstr)=i
579 ! deallocate(secstruc)
582 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
583 ! deallocate(secstruc)
586 write(iout,'(A20)')'Error reading FTORS'
587 ! deallocate(secstruc)
589 end subroutine secstrp2dihc
590 !-----------------------------------------------------------------------------
591 subroutine read_secstr_pred(jin,jout,errors)
593 ! implicit real*8 (a-h,o-z)
594 ! INCLUDE 'DIMENSIONS'
595 ! include 'COMMON.IOUNITS'
596 ! include 'COMMON.CHAIN'
597 !el character(len=1),dimension(nres) :: secstruc !(maxres)
598 !el COMMON/SECONDARYS/secstruc
600 character(len=80) :: line,line1 !,ucase
601 logical :: errflag,errors,blankline
604 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
607 read (jin,'(a)') line
608 write (jout,'(2a)') '> ',line(1:78)
610 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
611 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
612 ! end-groups in the input file "*.spred"
615 do while (index(line1,'$END').eq.0)
616 ! Override commented lines.
619 do while (.not.blankline)
621 call mykey(line,line1,ipos,blankline,errflag)
622 if (errflag) write (jout,'(2a)') &
623 'Error when reading sequence in line: ',line
624 errors=errors .or. errflag
625 if (.not. blankline .and. .not. errflag) then
628 !el if (iseq.le.maxres) then
629 if (line1(1:1).eq.'-' ) then
630 secstruc(iseq)=line1(1:1)
631 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
632 ( ucase(line1(1:1)).eq.'H' ) ) then
633 secstruc(iseq)=ucase(line1(1:1))
636 write (jout,1010) line1(1:1), iseq
641 !el write (jout,1000) iseq,maxres
644 do while (ipos1.le.iend)
649 !el if (iseq.le.maxres) then
650 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
651 secstruc(iseq)=line1(ipos1-1:ipos1-1)
652 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
653 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
654 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
657 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
662 !el write (jout,1000) iseq,maxres
669 read (jin,'(a)') line
670 write (jout,'(2a)') '> ',line(1:78)
674 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
676 !d check whether the found length of the chain is correct.
677 length_of_chain=iseq-1
678 if (length_of_chain .ne. nres) then
680 write (jout,'(a,i4,a,i4,a)') &
681 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
682 ,length_of_chain,') does not match with the number of residues (' &
687 1000 format('Error - the number of residues (',i4,&
688 ') has exceeded maximum (',i4,').')
689 1010 format ('Error - unrecognized secondary structure label',a4,&
692 end subroutine read_secstr_pred
694 !-----------------------------------------------------------------------------
696 !-----------------------------------------------------------------------------
701 use control_data, only:maxterm !,maxtor
705 use control, only: getenv_loc
707 ! Read the parameters of the probability distributions of the virtual-bond
708 ! valence angles and the side chains and energy parameters.
710 ! Important! Energy-term weights ARE NOT read here; they are read from the
711 ! main input file instead, because NO defaults have yet been set for these
714 ! implicit real*8 (a-h,o-z)
715 ! include 'DIMENSIONS'
720 ! include 'COMMON.IOUNITS'
721 ! include 'COMMON.CHAIN'
722 ! include 'COMMON.INTERACT'
723 ! include 'COMMON.GEO'
724 ! include 'COMMON.LOCAL'
725 ! include 'COMMON.TORSION'
726 ! include 'COMMON.SCCOR'
727 ! include 'COMMON.SCROT'
728 ! include 'COMMON.FFIELD'
729 ! include 'COMMON.NAMES'
730 ! include 'COMMON.SBRIDGE'
731 ! include 'COMMON.MD'
732 ! include 'COMMON.SETUP'
733 character(len=1) :: t1,t2,t3
734 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
735 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
736 logical :: lprint,LaTeX
737 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
738 real(kind=8),dimension(13) :: buse
739 character(len=3) :: lancuch !,ucase
741 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
742 integer :: maxinter,junk,kk,ii,ncatprotparm
743 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
744 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
745 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
746 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
748 integer :: ichir1,ichir2
749 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
750 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
751 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
753 ! For printing parameters after they are read set the following in the UNRES
756 ! setenv PRINT_PARM YES
758 ! To print parameters in LaTeX format rather than as ASCII tables:
762 call getenv_loc("PRINT_PARM",lancuch)
763 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
764 call getenv_loc("LATEX",lancuch)
765 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 dwa16=2.0d0**(1.0d0/6.0d0)
769 ! Assign virtual-bond length
772 vblinv2=vblinv*vblinv
774 ! Read the virtual-bond parameters, masses, and moments of inertia
775 ! and Stokes' radii of the peptide group and side chains
777 allocate(dsc(ntyp1)) !(ntyp1)
778 allocate(dsc_inv(ntyp1)) !(ntyp1)
779 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
780 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
781 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
782 allocate(nbondterm(ntyp)) !(ntyp)
783 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
784 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
785 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
786 allocate(long_r_sidechain(ntyp))
787 allocate(short_r_sidechain(ntyp))
792 allocate(msc(ntyp+1)) !(ntyp+1)
793 allocate(isc(ntyp+1)) !(ntyp+1)
794 allocate(restok(ntyp+1)) !(ntyp+1)
796 read (ibond,*) vbldp0,akp,mp,ip,pstok
799 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
800 dsc(i) = vbldsc0(1,i)
804 dsc_inv(i)=1.0D0/dsc(i)
808 allocate(msc(ntyp+1,5)) !(ntyp+1)
809 allocate(isc(ntyp+1,5)) !(ntyp+1)
810 allocate(restok(ntyp+1,5)) !(ntyp+1)
812 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
814 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
815 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
816 dsc(i) = vbldsc0(1,i)
820 dsc_inv(i)=1.0D0/dsc(i)
824 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
827 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
828 ! dsc(i) = vbldsc0_nucl(1,i)
832 ! dsc_inv(i)=1.0D0/dsc(i)
835 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
836 ! do i=1,ntyp_molec(2)
837 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
838 ! aksc_nucl(j,i),abond0_nucl(j,i),&
839 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
840 ! dsc(i) = vbldsc0(1,i)
844 ! dsc_inv(i)=1.0D0/dsc(i)
849 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
850 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
852 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
854 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
855 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
857 write (iout,'(13x,3f10.5)') &
858 vbldsc0(j,i),aksc(j,i),abond0(j,i)
863 read(iion,*) msc(i,5),restok(i,5)
864 print *,msc(i,5),restok(i,5)
868 read (iion,*) ncatprotparm
869 allocate(catprm(ncatprotparm,4))
871 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
874 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
875 !----------------------------------------------------
876 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
877 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
878 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
879 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
880 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
881 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
891 allocate(liptranene(ntyp))
892 !C reading lipid parameters
893 write (iout,*) "iliptranpar",iliptranpar
895 read(iliptranpar,*) pepliptran
898 read(iliptranpar,*) liptranene(i)
899 print *,liptranene(i)
905 ! Read the parameters of the probability distribution/energy expression
906 ! of the virtual-bond valence angles theta
909 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
910 (bthet(j,i,1,1),j=1,2)
911 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
912 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
913 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
917 athet(1,i,1,-1)=athet(1,i,1,1)
918 athet(2,i,1,-1)=athet(2,i,1,1)
919 bthet(1,i,1,-1)=-bthet(1,i,1,1)
920 bthet(2,i,1,-1)=-bthet(2,i,1,1)
921 athet(1,i,-1,1)=-athet(1,i,1,1)
922 athet(2,i,-1,1)=-athet(2,i,1,1)
923 bthet(1,i,-1,1)=bthet(1,i,1,1)
924 bthet(2,i,-1,1)=bthet(2,i,1,1)
928 athet(1,i,-1,-1)=athet(1,-i,1,1)
929 athet(2,i,-1,-1)=-athet(2,-i,1,1)
930 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
931 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
932 athet(1,i,-1,1)=athet(1,-i,1,1)
933 athet(2,i,-1,1)=-athet(2,-i,1,1)
934 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
935 bthet(2,i,-1,1)=bthet(2,-i,1,1)
936 athet(1,i,1,-1)=-athet(1,-i,1,1)
937 athet(2,i,1,-1)=athet(2,-i,1,1)
938 bthet(1,i,1,-1)=bthet(1,-i,1,1)
939 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
944 polthet(j,i)=polthet(j,-i)
947 gthet(j,i)=gthet(j,-i)
955 'Parameters of the virtual-bond valence angles:'
956 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
957 ' ATHETA0 ',' A1 ',' A2 ',&
960 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
961 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
963 write (iout,'(/a/9x,5a/79(1h-))') &
964 'Parameters of the expression for sigma(theta_c):',&
965 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
966 ' ALPH3 ',' SIGMA0C '
968 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
969 (polthet(j,i),j=0,3),sigc0(i)
971 write (iout,'(/a/9x,5a/79(1h-))') &
972 'Parameters of the second gaussian:',&
973 ' THETA0 ',' SIGMA0 ',' G1 ',&
976 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
977 sig0(i),(gthet(j,i),j=1,3)
981 'Parameters of the virtual-bond valence angles:'
982 write (iout,'(/a/9x,5a/79(1h-))') &
983 'Coefficients of expansion',&
984 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
985 ' b1*10^1 ',' b2*10^1 '
987 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
988 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
989 (10*bthet(j,i,1,1),j=1,2)
991 write (iout,'(/a/9x,5a/79(1h-))') &
992 'Parameters of the expression for sigma(theta_c):',&
993 ' alpha0 ',' alph1 ',' alph2 ',&
994 ' alhp3 ',' sigma0c '
996 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
997 (polthet(j,i),j=0,3),sigc0(i)
999 write (iout,'(/a/9x,5a/79(1h-))') &
1000 'Parameters of the second gaussian:',&
1001 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1004 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1005 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1011 ! Read the parameters of Utheta determined from ab initio surfaces
1012 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1014 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1015 ntheterm3,nsingle,ndouble
1016 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1018 !----------------------------------------------------
1019 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1020 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1021 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1022 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1023 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1024 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1025 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1026 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1027 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1028 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1029 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1030 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1031 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1032 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1033 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1034 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1035 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1036 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1037 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1038 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1039 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1040 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1041 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1042 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1044 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1046 ithetyp(i)=-ithetyp(-i)
1049 aa0thet(:,:,:,:)=0.0d0
1050 aathet(:,:,:,:,:)=0.0d0
1051 bbthet(:,:,:,:,:,:)=0.0d0
1052 ccthet(:,:,:,:,:,:)=0.0d0
1053 ddthet(:,:,:,:,:,:)=0.0d0
1054 eethet(:,:,:,:,:,:)=0.0d0
1055 ffthet(:,:,:,:,:,:,:)=0.0d0
1056 ggthet(:,:,:,:,:,:,:)=0.0d0
1058 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1060 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1061 ! VAR:1=non-glicyne non-proline 2=proline
1062 ! VAR:negative values for D-aminoacid
1064 do j=-nthetyp,nthetyp
1065 do k=-nthetyp,nthetyp
1066 read (ithep,'(6a)',end=111,err=111) res1
1067 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1068 ! VAR: aa0thet is variable describing the average value of Foureir
1069 ! VAR: expansion series
1070 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1071 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1072 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1073 read (ithep,*,end=111,err=111) &
1074 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1075 read (ithep,*,end=111,err=111) &
1076 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1077 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1078 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1079 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1081 read (ithep,*,end=111,err=111) &
1082 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1083 ffthet(lll,llll,ll,i,j,k,iblock),&
1084 ggthet(llll,lll,ll,i,j,k,iblock),&
1085 ggthet(lll,llll,ll,i,j,k,iblock),&
1086 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1091 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1092 ! coefficients of theta-and-gamma-dependent terms are zero.
1093 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1094 ! RECOMENTDED AFTER VERSION 3.3)
1098 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1099 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1101 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1102 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1105 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1107 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1110 ! AND COMMENT THE LOOPS BELOW
1114 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1115 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1117 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1118 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1121 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1123 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1128 ! Substitution for D aminoacids from symmetry.
1131 do j=-nthetyp,nthetyp
1132 do k=-nthetyp,nthetyp
1133 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1135 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1139 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1140 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1141 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1142 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1148 ffthet(llll,lll,ll,i,j,k,iblock)= &
1149 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1150 ffthet(lll,llll,ll,i,j,k,iblock)= &
1151 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1152 ggthet(llll,lll,ll,i,j,k,iblock)= &
1153 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1154 ggthet(lll,llll,ll,i,j,k,iblock)= &
1155 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1164 ! Control printout of the coefficients of virtual-bond-angle potentials
1167 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1172 write (iout,'(//4a)') &
1173 'Type ',onelett(i),onelett(j),onelett(k)
1174 write (iout,'(//a,10x,a)') " l","a[l]"
1175 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1176 write (iout,'(i2,1pe15.5)') &
1177 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1179 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1180 "b",l,"c",l,"d",l,"e",l
1182 write (iout,'(i2,4(1pe15.5))') m,&
1183 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1184 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1188 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1189 "f+",l,"f-",l,"g+",l,"g-",l
1192 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1193 ffthet(n,m,l,i,j,k,iblock),&
1194 ffthet(m,n,l,i,j,k,iblock),&
1195 ggthet(n,m,l,i,j,k,iblock),&
1196 ggthet(m,n,l,i,j,k,iblock)
1206 write (2,*) "Start reading THETA_PDB",ithep_pdb
1208 ! write (2,*) 'i=',i
1209 read (ithep_pdb,*,err=111,end=111) &
1210 a0thet(i),(athet(j,i,1,1),j=1,2),&
1211 (bthet(j,i,1,1),j=1,2)
1212 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1213 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1214 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1215 sigc0(i)=sigc0(i)**2
1218 athet(1,i,1,-1)=athet(1,i,1,1)
1219 athet(2,i,1,-1)=athet(2,i,1,1)
1220 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1221 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1222 athet(1,i,-1,1)=-athet(1,i,1,1)
1223 athet(2,i,-1,1)=-athet(2,i,1,1)
1224 bthet(1,i,-1,1)=bthet(1,i,1,1)
1225 bthet(2,i,-1,1)=bthet(2,i,1,1)
1228 a0thet(i)=a0thet(-i)
1229 athet(1,i,-1,-1)=athet(1,-i,1,1)
1230 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1231 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1232 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1233 athet(1,i,-1,1)=athet(1,-i,1,1)
1234 athet(2,i,-1,1)=-athet(2,-i,1,1)
1235 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1236 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1237 athet(1,i,1,-1)=-athet(1,-i,1,1)
1238 athet(2,i,1,-1)=athet(2,-i,1,1)
1239 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1240 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1241 theta0(i)=theta0(-i)
1245 polthet(j,i)=polthet(j,-i)
1248 gthet(j,i)=gthet(j,-i)
1251 write (2,*) "End reading THETA_PDB"
1255 !--------------- Reading theta parameters for nucleic acid-------
1256 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1257 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1258 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1259 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1260 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1261 nthetyp_nucl+1,nthetyp_nucl+1))
1262 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1263 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1264 nthetyp_nucl+1,nthetyp_nucl+1))
1265 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1266 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1267 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1268 nthetyp_nucl+1,nthetyp_nucl+1))
1269 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1270 nthetyp_nucl+1,nthetyp_nucl+1))
1271 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1272 nthetyp_nucl+1,nthetyp_nucl+1))
1273 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1274 nthetyp_nucl+1,nthetyp_nucl+1))
1275 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1276 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1277 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1278 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1279 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1280 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1282 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1283 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1285 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1287 aa0thet_nucl(:,:,:)=0.0d0
1288 aathet_nucl(:,:,:,:)=0.0d0
1289 bbthet_nucl(:,:,:,:,:)=0.0d0
1290 ccthet_nucl(:,:,:,:,:)=0.0d0
1291 ddthet_nucl(:,:,:,:,:)=0.0d0
1292 eethet_nucl(:,:,:,:,:)=0.0d0
1293 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1294 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1299 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1300 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1301 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1302 read (ithep_nucl,*,end=111,err=111) &
1303 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1304 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1305 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1306 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1307 read (ithep_nucl,*,end=111,err=111) &
1308 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1309 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1310 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1315 !-------------------------------------------
1316 allocate(nlob(ntyp1)) !(ntyp1)
1317 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1318 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1319 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1326 gaussc(:,:,:,:)=0.0D0
1330 ! Read the parameters of the probability distribution/energy expression
1331 ! of the side chains.
1334 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1338 dsc_inv(i)=1.0D0/dsc(i)
1349 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1350 ((blower(k,l,1),l=1,k),k=1,3)
1351 censc(1,1,-i)=censc(1,1,i)
1352 censc(2,1,-i)=censc(2,1,i)
1353 censc(3,1,-i)=-censc(3,1,i)
1355 read (irotam,*,end=112,err=112) bsc(j,i)
1356 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1357 ((blower(k,l,j),l=1,k),k=1,3)
1358 censc(1,j,-i)=censc(1,j,i)
1359 censc(2,j,-i)=censc(2,j,i)
1360 censc(3,j,-i)=-censc(3,j,i)
1361 ! BSC is amplitude of Gaussian
1368 akl=akl+blower(k,m,j)*blower(l,m,j)
1372 if (((k.eq.3).and.(l.ne.3)) &
1373 .or.((l.eq.3).and.(k.ne.3))) then
1374 gaussc(k,l,j,-i)=-akl
1375 gaussc(l,k,j,-i)=-akl
1377 gaussc(k,l,j,-i)=akl
1378 gaussc(l,k,j,-i)=akl
1387 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1390 if (nlobi.gt.0) then
1392 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1393 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1394 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1395 'log h',(bsc(j,i),j=1,nlobi)
1396 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1397 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1399 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1400 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1403 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1404 write (iout,'(a,f10.4,4(16x,f10.4))') &
1405 'Center ',(bsc(j,i),j=1,nlobi)
1406 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1415 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1416 ! added by Urszula Kozlowska 07/11/2007
1418 !el Maximum number of SC local term fitting function coefficiants
1419 !el integer,parameter :: maxsccoef=65
1421 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1424 read (irotam,*,end=112,err=112)
1426 read (irotam,*,end=112,err=112)
1429 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1433 !---------reading nucleic acid parameters for rotamers-------------------
1434 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1435 do i=1,ntyp_molec(2)
1436 read (irotam_nucl,*,end=112,err=112)
1438 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1444 write (iout,*) "Base rotamer parameters"
1445 do i=1,ntyp_molec(2)
1446 write (iout,'(a)') restyp(i,2)
1447 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1452 ! Read the parameters of the probability distribution/energy expression
1453 ! of the side chains.
1455 write (2,*) "Start reading ROTAM_PDB"
1457 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1461 dsc_inv(i)=1.0D0/dsc(i)
1472 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1473 ((blower(k,l,1),l=1,k),k=1,3)
1475 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1476 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1477 ((blower(k,l,j),l=1,k),k=1,3)
1484 akl=akl+blower(k,m,j)*blower(l,m,j)
1494 write (2,*) "End reading ROTAM_PDB"
1500 ! Read torsional parameters in old format
1502 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1504 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1505 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1506 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1508 !el from energy module--------
1509 allocate(v1(nterm_old,ntortyp,ntortyp))
1510 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1511 !el---------------------------
1516 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1522 write (iout,'(/a/)') 'Torsional constants:'
1525 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1526 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1532 ! Read torsional parameters
1534 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1535 read (itorp,*,end=113,err=113) ntortyp
1536 !el from energy module---------
1537 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1538 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1540 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1541 allocate(vlor2(maxlor,ntortyp,ntortyp))
1542 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1543 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1545 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1546 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1547 !el---------------------------
1550 !el---------------------------
1552 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1554 itortyp(i)=-itortyp(-i)
1556 itortyp(ntyp1)=ntortyp
1557 itortyp(-ntyp1)=-ntortyp
1559 write (iout,*) 'ntortyp',ntortyp
1561 do j=-ntortyp+1,ntortyp-1
1562 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1564 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1565 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1568 do k=1,nterm(i,j,iblock)
1569 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1571 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1572 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1573 v0ij=v0ij+si*v1(k,i,j,iblock)
1575 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1576 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1577 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1579 do k=1,nlor(i,j,iblock)
1580 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1581 vlor2(k,i,j),vlor3(k,i,j)
1582 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1585 v0(-i,-j,iblock)=v0ij
1591 write (iout,'(/a/)') 'Torsional constants:'
1593 do i=-ntortyp,ntortyp
1594 do j=-ntortyp,ntortyp
1595 write (iout,*) 'ityp',i,' jtyp',j
1596 write (iout,*) 'Fourier constants'
1597 do k=1,nterm(i,j,iblock)
1598 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1601 write (iout,*) 'Lorenz constants'
1602 do k=1,nlor(i,j,iblock)
1603 write (iout,'(3(1pe15.5))') &
1604 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1610 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1612 ! 6/23/01 Read parameters for double torsionals
1614 !el from energy module------------
1615 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1616 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1617 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1618 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1619 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1620 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1621 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1622 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1623 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1624 !---------------------------------
1628 do j=-ntortyp+1,ntortyp-1
1629 do k=-ntortyp+1,ntortyp-1
1630 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1631 ! write (iout,*) "OK onelett",
1634 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1635 .or. t3.ne.toronelet(k)) then
1636 write (iout,*) "Error in double torsional parameter file",&
1639 call MPI_Finalize(Ierror)
1641 stop "Error in double torsional parameter file"
1643 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1644 ntermd_2(i,j,k,iblock)
1645 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1646 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1647 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1648 ntermd_1(i,j,k,iblock))
1649 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1650 ntermd_1(i,j,k,iblock))
1651 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1652 ntermd_1(i,j,k,iblock))
1653 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1654 ntermd_1(i,j,k,iblock))
1655 ! Martix of D parameters for one dimesional foureir series
1656 do l=1,ntermd_1(i,j,k,iblock)
1657 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1658 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1659 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1660 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1661 ! write(iout,*) "whcodze" ,
1662 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1664 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1665 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1666 v2s(m,l,i,j,k,iblock),&
1667 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1668 ! Martix of D parameters for two dimesional fourier series
1669 do l=1,ntermd_2(i,j,k,iblock)
1671 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1672 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1673 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1674 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1683 write (iout,*) 'Constants for double torsionals'
1686 do j=-ntortyp+1,ntortyp-1
1687 do k=-ntortyp+1,ntortyp-1
1688 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1689 ' nsingle',ntermd_1(i,j,k,iblock),&
1690 ' ndouble',ntermd_2(i,j,k,iblock)
1692 write (iout,*) 'Single angles:'
1693 do l=1,ntermd_1(i,j,k,iblock)
1694 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1695 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1696 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1697 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1700 write (iout,*) 'Pairs of angles:'
1701 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1702 do l=1,ntermd_2(i,j,k,iblock)
1703 write (iout,'(i5,20f10.5)') &
1704 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1707 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1708 do l=1,ntermd_2(i,j,k,iblock)
1709 write (iout,'(i5,20f10.5)') &
1710 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1711 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1720 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1721 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
1722 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
1723 !el from energy module---------
1724 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1725 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1727 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1728 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
1729 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
1730 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1732 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
1733 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1734 !el---------------------------
1737 !el--------------------
1738 read (itorp_nucl,*,end=113,err=113) &
1739 (itortyp_nucl(i),i=1,ntyp_molec(2))
1740 ! print *,itortyp_nucl(:)
1741 !c write (iout,*) 'ntortyp',ntortyp
1744 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
1745 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
1748 do k=1,nterm_nucl(i,j)
1749 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
1750 v0ij=v0ij+si*v1_nucl(k,i,j)
1753 do k=1,nlor_nucl(i,j)
1754 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
1755 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
1756 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
1762 ! Read of Side-chain backbone correlation parameters
1763 ! Modified 11 May 2012 by Adasko
1766 read (isccor,*,end=119,err=119) nsccortyp
1768 !el from module energy-------------
1769 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1770 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1771 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1772 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1773 !-----------------------------------
1775 !el from module energy-------------
1776 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1778 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1780 isccortyp(i)=-isccortyp(-i)
1782 iscprol=isccortyp(20)
1783 ! write (iout,*) 'ntortyp',ntortyp
1785 !c maxinter is maximum interaction sites
1786 !el from module energy---------
1787 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1788 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1789 -nsccortyp:nsccortyp))
1790 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1791 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1792 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1793 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1794 !-----------------------------------
1798 read (isccor,*,end=119,err=119) &
1799 nterm_sccor(i,j),nlor_sccor(i,j)
1805 nterm_sccor(-i,j)=nterm_sccor(i,j)
1806 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1807 nterm_sccor(i,-j)=nterm_sccor(i,j)
1808 do k=1,nterm_sccor(i,j)
1809 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1811 if (j.eq.iscprol) then
1812 if (i.eq.isccortyp(10)) then
1813 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1814 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1816 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1817 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1818 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1819 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1820 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1821 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1822 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1823 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1826 if (i.eq.isccortyp(10)) then
1827 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1828 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1830 if (j.eq.isccortyp(10)) then
1831 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1832 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1834 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1835 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1836 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1837 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1838 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1839 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1843 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1844 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1845 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1846 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1849 do k=1,nlor_sccor(i,j)
1850 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1851 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1852 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1853 (1+vlor3sccor(k,i,j)**2)
1855 v0sccor(l,i,j)=v0ijsccor
1856 v0sccor(l,-i,j)=v0ijsccor1
1857 v0sccor(l,i,-j)=v0ijsccor2
1858 v0sccor(l,-i,-j)=v0ijsccor3
1864 !el from module energy-------------
1865 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1867 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1868 ! write (iout,*) 'ntortyp',ntortyp
1870 !c maxinter is maximum interaction sites
1871 !el from module energy---------
1872 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1873 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1874 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1875 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1876 !-----------------------------------
1880 read (isccor,*,end=119,err=119) &
1881 nterm_sccor(i,j),nlor_sccor(i,j)
1885 do k=1,nterm_sccor(i,j)
1886 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1888 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1891 do k=1,nlor_sccor(i,j)
1892 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1893 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1894 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1895 (1+vlor3sccor(k,i,j)**2)
1897 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1906 write (iout,'(/a/)') 'Torsional constants:'
1909 write (iout,*) 'ityp',i,' jtyp',j
1910 write (iout,*) 'Fourier constants'
1911 do k=1,nterm_sccor(i,j)
1912 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1914 write (iout,*) 'Lorenz constants'
1915 do k=1,nlor_sccor(i,j)
1916 write (iout,'(3(1pe15.5))') &
1917 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1924 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1925 ! interaction energy of the Gly, Ala, and Pro prototypes.
1930 write (iout,*) "Coefficients of the cumulants"
1932 read (ifourier,*) nloctyp
1933 !write(iout,*) "nloctyp",nloctyp
1934 !el from module energy-------
1935 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1936 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1937 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1938 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1939 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1940 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1941 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1942 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1946 !--------------------------------
1949 read (ifourier,*,end=115,err=115)
1950 read (ifourier,*,end=115,err=115) (buse(ii),ii=1,13)
1952 write (iout,*) 'Type',i
1953 write (iout,'(a,i2,a,f10.5)') ('buse(',ii,')=',buse(ii),ii=1,13)
1961 B1tilde(1,i) = buse(3)
1962 B1tilde(2,i) =-buse(5)
1963 B1tilde(1,-i) =-buse(3)
1964 B1tilde(2,-i) =buse(5)
1965 ! buse1tilde(1,i)=0.0d0
1966 ! buse1tilde(2,i)=0.0d0
1986 Ctilde(1,1,i)=buse(7)
1987 Ctilde(1,2,i)=buse(9)
1988 Ctilde(2,1,i)=-buse(9)
1989 Ctilde(2,2,i)=buse(7)
1990 Ctilde(1,1,-i)=buse(7)
1991 Ctilde(1,2,-i)=-buse(9)
1992 Ctilde(2,1,-i)=buse(9)
1993 Ctilde(2,2,-i)=buse(7)
1995 ! Ctilde(1,1,i)=0.0d0
1996 ! Ctilde(1,2,i)=0.0d0
1997 ! Ctilde(2,1,i)=0.0d0
1998 ! Ctilde(2,2,i)=0.0d0
2011 Dtilde(1,1,i)=buse(6)
2012 Dtilde(1,2,i)=buse(8)
2013 Dtilde(2,1,i)=-buse(8)
2014 Dtilde(2,2,i)=buse(6)
2015 Dtilde(1,1,-i)=buse(6)
2016 Dtilde(1,2,-i)=-buse(8)
2017 Dtilde(2,1,-i)=buse(8)
2018 Dtilde(2,2,-i)=buse(6)
2020 ! Dtilde(1,1,i)=0.0d0
2021 ! Dtilde(1,2,i)=0.0d0
2022 ! Dtilde(2,1,i)=0.0d0
2023 ! Dtilde(2,2,i)=0.0d0
2024 EE(1,1,i)= buse(10)+buse(11)
2025 EE(2,2,i)=-buse(10)+buse(11)
2026 EE(2,1,i)= buse(12)-buse(13)
2027 EE(1,2,i)= buse(12)+buse(13)
2028 EE(1,1,-i)= buse(10)+buse(11)
2029 EE(2,2,-i)=-buse(10)+buse(11)
2030 EE(2,1,-i)=-buse(12)+buse(13)
2031 EE(1,2,-i)=-buse(12)-buse(13)
2037 ! ee(2,1,i)=ee(1,2,i)
2041 write (iout,*) 'Type',i
2043 write(iout,*) B1(1,i),B1(2,i)
2045 write(iout,*) B2(1,i),B2(2,i)
2048 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
2052 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
2056 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
2061 ! Read electrostatic-interaction parameters
2066 write (iout,'(/a)') 'Electrostatic interaction constants:'
2067 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2068 'IT','JT','APP','BPP','AEL6','AEL3'
2070 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2071 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2072 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2073 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2078 app (i,j)=epp(i,j)*rri*rri
2079 bpp (i,j)=-2.0D0*epp(i,j)*rri
2080 ael6(i,j)=elpp6(i,j)*4.2D0**6
2081 ael3(i,j)=elpp3(i,j)*4.2D0**3
2083 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2089 ! Read side-chain interaction parameters.
2091 !el from module energy - COMMON.INTERACT-------
2092 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2093 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2094 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2096 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2097 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2098 allocate(epslip(ntyp,ntyp))
2106 !--------------------------------
2108 read (isidep,*,end=117,err=117) ipot,expon
2109 if (ipot.lt.1 .or. ipot.gt.5) then
2110 write (iout,'(2a)') 'Error while reading SC interaction',&
2111 'potential file - unknown potential type.'
2113 call MPI_Finalize(Ierror)
2119 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2120 ', exponents are ',expon,2*expon
2121 ! goto (10,20,30,30,40) ipot
2123 !----------------------- LJ potential ---------------------------------
2125 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2126 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2127 (sigma0(i),i=1,ntyp)
2129 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2130 write (iout,'(a/)') 'The epsilon array:'
2131 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2132 write (iout,'(/a)') 'One-body parameters:'
2133 write (iout,'(a,4x,a)') 'residue','sigma'
2134 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2137 !----------------------- LJK potential --------------------------------
2139 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2140 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2141 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2143 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2144 write (iout,'(a/)') 'The epsilon array:'
2145 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2146 write (iout,'(/a)') 'One-body parameters:'
2147 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2148 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2152 !---------------------- GB or BP potential -----------------------------
2156 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2158 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2159 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2160 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2161 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2163 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2166 ! For the GB potential convert sigma'**2 into chi'
2169 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2173 write (iout,'(/a/)') 'Parameters of the BP potential:'
2174 write (iout,'(a/)') 'The epsilon array:'
2175 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2176 write (iout,'(/a)') 'One-body parameters:'
2177 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2179 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2180 chip(i),alp(i),i=1,ntyp)
2183 !--------------------- GBV potential -----------------------------------
2185 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2186 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2187 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2188 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2190 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2191 write (iout,'(a/)') 'The epsilon array:'
2192 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2193 write (iout,'(/a)') 'One-body parameters:'
2194 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2195 's||/s_|_^2',' chip ',' alph '
2196 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2197 sigii(i),chip(i),alp(i),i=1,ntyp)
2200 write(iout,*)"Wrong ipot"
2206 !-----------------------------------------------------------------------
2207 ! Calculate the "working" parameters of SC interactions.
2209 !el from module energy - COMMON.INTERACT-------
2210 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2211 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2212 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2213 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2215 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2228 sc_aa_tube_par(:)=0.0d0
2229 sc_bb_tube_par(:)=0.0d0
2231 !--------------------------------
2236 epslip(i,j)=epslip(j,i)
2241 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2242 sigma(j,i)=sigma(i,j)
2243 rs0(i,j)=dwa16*sigma(i,j)
2247 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2248 'Working parameters of the SC interactions:',&
2249 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2254 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2263 sigeps=dsign(1.0D0,epsij)
2265 aa_aq(i,j)=epsij*rrij*rrij
2266 bb_aq(i,j)=-sigeps*epsij*rrij
2267 aa_aq(j,i)=aa_aq(i,j)
2268 bb_aq(j,i)=bb_aq(i,j)
2269 epsijlip=epslip(i,j)
2270 sigeps=dsign(1.0D0,epsijlip)
2271 epsijlip=dabs(epsijlip)
2272 aa_lip(i,j)=epsijlip*rrij*rrij
2273 bb_lip(i,j)=-sigeps*epsijlip*rrij
2274 aa_lip(j,i)=aa_lip(i,j)
2275 bb_lip(j,i)=bb_lip(i,j)
2276 !C write(iout,*) aa_lip
2278 sigt1sq=sigma0(i)**2
2279 sigt2sq=sigma0(j)**2
2282 ratsig1=sigt2sq/sigt1sq
2283 ratsig2=1.0D0/ratsig1
2284 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2285 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2286 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2290 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2291 sigmaii(i,j)=rsum_max
2292 sigmaii(j,i)=rsum_max
2294 ! sigmaii(i,j)=r0(i,j)
2295 ! sigmaii(j,i)=r0(i,j)
2297 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2298 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2299 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2300 augm(i,j)=epsij*r_augm**(2*expon)
2301 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2308 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2309 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2310 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2315 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2316 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2317 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2318 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2319 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2320 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2321 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2322 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2323 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2324 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2325 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2326 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2327 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2328 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2329 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2330 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2339 read (isidep_nucl,*) ipot_nucl
2340 ! print *,"TU?!",ipot_nucl
2341 if (ipot_nucl.eq.1) then
2342 do i=1,ntyp_molec(2)
2343 do j=i,ntyp_molec(2)
2344 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2345 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2349 do i=1,ntyp_molec(2)
2350 do j=i,ntyp_molec(2)
2351 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2352 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2353 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2357 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2358 do i=1,ntyp_molec(2)
2359 do j=i,ntyp_molec(2)
2360 rrij=sigma_nucl(i,j)
2364 epsij=4*eps_nucl(i,j)
2365 sigeps=dsign(1.0D0,epsij)
2367 aa_nucl(i,j)=epsij*rrij*rrij
2368 bb_nucl(i,j)=-sigeps*epsij*rrij
2369 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2370 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2371 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2372 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2373 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2374 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2377 aa_nucl(i,j)=aa_nucl(j,i)
2378 bb_nucl(i,j)=bb_nucl(j,i)
2379 ael3_nucl(i,j)=ael3_nucl(j,i)
2380 ael6_nucl(i,j)=ael6_nucl(j,i)
2381 ael63_nucl(i,j)=ael63_nucl(j,i)
2382 ael32_nucl(i,j)=ael32_nucl(j,i)
2383 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2384 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2385 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2386 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2387 eps_nucl(i,j)=eps_nucl(j,i)
2388 sigma_nucl(i,j)=sigma_nucl(j,i)
2389 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2393 write(iout,*) "tube param"
2394 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2395 ccavtubpep,dcavtubpep,tubetranenepep
2396 sigmapeptube=sigmapeptube**6
2397 sigeps=dsign(1.0D0,epspeptube)
2398 epspeptube=dabs(epspeptube)
2399 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2400 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2401 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2403 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2404 ccavtub(i),dcavtub(i),tubetranene(i)
2405 sigmasctube=sigmasctube**6
2406 sigeps=dsign(1.0D0,epssctube)
2407 epssctube=dabs(epssctube)
2408 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2409 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2410 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2412 !-----------------READING SC BASE POTENTIALS-----------------------------
2413 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2414 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2415 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2416 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2417 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2418 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2419 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2420 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2421 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2422 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2423 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2424 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2425 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2426 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2427 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2428 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2431 do i=1,ntyp_molec(1)
2432 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2433 write (*,*) "Im in ", i, " ", j
2434 read(isidep_scbase,*) &
2435 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2436 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2437 write(*,*) "eps",eps_scbase(i,j)
2438 read(isidep_scbase,*) &
2439 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2440 chis_scbase(i,j,1),chis_scbase(i,j,2)
2441 read(isidep_scbase,*) &
2442 dhead_scbasei(i,j), &
2443 dhead_scbasej(i,j), &
2444 rborn_scbasei(i,j),rborn_scbasej(i,j)
2445 read(isidep_scbase,*) &
2446 (wdipdip_scbase(k,i,j),k=1,3), &
2447 (wqdip_scbase(k,i,j),k=1,2)
2448 read(isidep_scbase,*) &
2449 alphapol_scbase(i,j), &
2450 epsintab_scbase(i,j)
2453 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2454 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2456 do i=1,ntyp_molec(1)
2457 do j=1,ntyp_molec(2)-1
2458 epsij=eps_scbase(i,j)
2459 rrij=sigma_scbase(i,j)
2464 sigeps=dsign(1.0D0,epsij)
2466 aa_scbase(i,j)=epsij*rrij*rrij
2467 bb_scbase(i,j)=-sigeps*epsij*rrij
2470 !-----------------READING PEP BASE POTENTIALS-------------------
2471 allocate(eps_pepbase(ntyp_molec(2)))
2472 allocate(sigma_pepbase(ntyp_molec(2)))
2473 allocate(chi_pepbase(ntyp_molec(2),2))
2474 allocate(chipp_pepbase(ntyp_molec(2),2))
2475 allocate(alphasur_pepbase(4,ntyp_molec(2)))
2476 allocate(sigmap1_pepbase(ntyp_molec(2)))
2477 allocate(sigmap2_pepbase(ntyp_molec(2)))
2478 allocate(chis_pepbase(ntyp_molec(2),2))
2479 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
2482 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2483 write (*,*) "Im in ", i, " ", j
2484 read(isidep_pepbase,*) &
2485 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
2486 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
2487 write(*,*) "eps",eps_pepbase(j)
2488 read(isidep_pepbase,*) &
2489 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
2490 chis_pepbase(j,1),chis_pepbase(j,2)
2491 read(isidep_pepbase,*) &
2492 (wdipdip_pepbase(k,j),k=1,3)
2494 allocate(aa_pepbase(ntyp_molec(2)))
2495 allocate(bb_pepbase(ntyp_molec(2)))
2497 do j=1,ntyp_molec(2)-1
2498 epsij=eps_pepbase(j)
2499 rrij=sigma_pepbase(j)
2504 sigeps=dsign(1.0D0,epsij)
2506 aa_pepbase(j)=epsij*rrij*rrij
2507 bb_pepbase(j)=-sigeps*epsij*rrij
2509 !--------------READING SC PHOSPHATE-------------------------------------
2510 allocate(eps_scpho(ntyp_molec(1)))
2511 allocate(sigma_scpho(ntyp_molec(1)))
2512 allocate(chi_scpho(ntyp_molec(1),2))
2513 allocate(chipp_scpho(ntyp_molec(1),2))
2514 allocate(alphasur_scpho(4,ntyp_molec(1)))
2515 allocate(sigmap1_scpho(ntyp_molec(1)))
2516 allocate(sigmap2_scpho(ntyp_molec(1)))
2517 allocate(chis_scpho(ntyp_molec(1),2))
2518 allocate(wqq_scpho(ntyp_molec(1)))
2519 allocate(wqdip_scpho(2,ntyp_molec(1)))
2520 allocate(alphapol_scpho(ntyp_molec(1)))
2521 allocate(epsintab_scpho(ntyp_molec(1)))
2522 allocate(dhead_scphoi(ntyp_molec(1)))
2523 allocate(rborn_scphoi(ntyp_molec(1)))
2524 allocate(rborn_scphoj(ntyp_molec(1)))
2525 allocate(alphi_scpho(ntyp_molec(1)))
2529 do j=1,ntyp_molec(1) ! without U then we will take T for U
2530 write (*,*) "Im in scpho ", i, " ", j
2531 read(isidep_scpho,*) &
2532 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
2533 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
2534 write(*,*) "eps",eps_scpho(j)
2535 read(isidep_scpho,*) &
2536 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
2537 chis_scpho(j,1),chis_scpho(j,2)
2538 read(isidep_scpho,*) &
2539 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
2540 read(isidep_scpho,*) &
2541 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
2545 allocate(aa_scpho(ntyp_molec(1)))
2546 allocate(bb_scpho(ntyp_molec(1)))
2548 do j=1,ntyp_molec(1)
2555 sigeps=dsign(1.0D0,epsij)
2557 aa_scpho(j)=epsij*rrij*rrij
2558 bb_scpho(j)=-sigeps*epsij*rrij
2562 read(isidep_peppho,*) &
2563 eps_peppho,sigma_peppho
2564 read(isidep_peppho,*) &
2565 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
2566 read(isidep_peppho,*) &
2567 (wqdip_peppho(k),k=1,2)
2575 sigeps=dsign(1.0D0,epsij)
2577 aa_peppho=epsij*rrij*rrij
2578 bb_peppho=-sigeps*epsij*rrij
2581 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2586 ! Define the SC-p interaction constants (hard-coded; old style)
2589 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2591 ! aad(i,1)=0.3D0*4.0D0**12
2592 ! Following line for constants currently implemented
2593 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2594 aad(i,1)=1.5D0*4.0D0**12
2595 ! aad(i,1)=0.17D0*5.6D0**12
2597 ! "Soft" SC-p repulsion
2599 ! Following line for constants currently implemented
2600 ! aad(i,1)=0.3D0*4.0D0**6
2601 ! "Hard" SC-p repulsion
2602 bad(i,1)=3.0D0*4.0D0**6
2603 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2612 ! 8/9/01 Read the SC-p interaction constants from file
2615 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2618 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2619 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2620 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2621 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2625 write (iout,*) "Parameters of SC-p interactions:"
2627 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2628 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2633 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2635 do i=1,ntyp_molec(2)
2636 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
2638 do i=1,ntyp_molec(2)
2639 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2640 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2642 r0pp=1.12246204830937298142*5.16158
2648 ! Define the constants of the disulfide bridge
2652 ! Old arbitrary potential - commented out.
2657 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2658 ! energy surface of diethyl disulfide.
2659 ! A. Liwo and U. Kozlowska, 11/24/03
2676 write (iout,'(/a)') "Disulfide bridge parameters:"
2677 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2678 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2679 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2680 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2684 111 write (iout,*) "Error reading bending energy parameters."
2686 112 write (iout,*) "Error reading rotamer energy parameters."
2688 113 write (iout,*) "Error reading torsional energy parameters."
2690 114 write (iout,*) "Error reading double torsional energy parameters."
2692 115 write (iout,*) &
2693 "Error reading cumulant (multibody energy) parameters."
2695 116 write (iout,*) "Error reading electrostatic energy parameters."
2697 117 write (iout,*) "Error reading side chain interaction parameters."
2699 118 write (iout,*) "Error reading SCp interaction parameters."
2701 119 write (iout,*) "Error reading SCCOR parameters"
2704 call MPI_Finalize(Ierror)
2708 end subroutine parmread
2710 !-----------------------------------------------------------------------------
2712 !-----------------------------------------------------------------------------
2713 subroutine printmat(ldim,m,n,iout,key,a)
2716 character(len=3),dimension(n) :: key
2717 real(kind=8),dimension(ldim,n) :: a
2719 integer :: i,j,k,m,iout,nlim
2723 write (iout,1000) (key(k),k=i,nlim)
2725 1000 format (/5x,8(6x,a3))
2726 1020 format (/80(1h-)/)
2728 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2731 1010 format (a3,2x,8(f9.4))
2733 end subroutine printmat
2734 !-----------------------------------------------------------------------------
2736 !-----------------------------------------------------------------------------
2738 ! Read the PDB file and convert the peptide geometry into virtual-chain
2741 use energy_data, only: itype,istype
2745 use control, only: rescode,sugarcode
2746 ! implicit real*8 (a-h,o-z)
2747 ! include 'DIMENSIONS'
2748 ! include 'COMMON.LOCAL'
2749 ! include 'COMMON.VAR'
2750 ! include 'COMMON.CHAIN'
2751 ! include 'COMMON.INTERACT'
2752 ! include 'COMMON.IOUNITS'
2753 ! include 'COMMON.GEO'
2754 ! include 'COMMON.NAMES'
2755 ! include 'COMMON.CONTROL'
2756 ! include 'COMMON.DISTFIT'
2757 ! include 'COMMON.SETUP'
2758 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
2760 logical :: lprn=.true.,fail
2761 real(kind=8),dimension(3) :: e1,e2,e3
2762 real(kind=8) :: dcj,efree_temp
2763 character(len=3) :: seq,res
2764 character(len=5) :: atom
2765 character(len=80) :: card
2766 real(kind=8),dimension(3,20) :: sccor
2767 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
2768 integer :: isugar,molecprev,firstion
2769 character*1 :: sugar
2771 real(kind=8),dimension(3) :: ccc
2773 integer,dimension(2,maxres/3) :: hfrag_alloc
2774 integer,dimension(4,maxres/3) :: bfrag_alloc
2775 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2776 real(kind=8),dimension(:,:), allocatable :: c_temporary
2777 integer,dimension(:,:) , allocatable :: itype_temporary
2778 integer,dimension(:),allocatable :: istype_temp
2785 ! write (2,*) "UNRES_PDB",unres_pdb
2805 !-----------------------------
2806 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2807 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2808 if(.not. allocated(istype)) allocate(istype(maxres))
2810 read (ipdbin,'(a80)',end=10) card
2811 ! write (iout,'(a)') card
2812 if (card(:5).eq.'HELIX') then
2815 read(card(22:25),*) hfrag(1,nhfrag)
2816 read(card(34:37),*) hfrag(2,nhfrag)
2818 if (card(:5).eq.'SHEET') then
2821 read(card(24:26),*) bfrag(1,nbfrag)
2822 read(card(35:37),*) bfrag(2,nbfrag)
2823 !rc----------------------------------------
2824 !rc to be corrected !!!
2825 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2826 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2827 !rc----------------------------------------
2829 if (card(:3).eq.'END') then
2831 else if (card(:3).eq.'TER') then
2836 itype(ires_old,molecule)=ntyp1_molec(molecule)
2837 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
2838 nres_molec(molecule)=nres_molec(molecule)+2
2840 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2843 dc(j,ires)=sccor(j,iii)
2846 call sccenter(ires,iii,sccor)
2852 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2853 ! Fish out the ATOM cards.
2854 ! write(iout,*) 'card',card(1:20)
2855 if (index(card(1:4),'ATOM').gt.0) then
2856 read (card(12:16),*) atom
2857 ! write (iout,*) "! ",atom," !",ires
2858 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2859 read (card(23:26),*) ires
2860 read (card(18:20),'(a3)') res
2861 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2862 ! & " ires_old",ires_old
2863 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2864 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2865 if (ires-ishift+ishift1.ne.ires_old) then
2866 ! Calculate the CM of the preceding residue.
2867 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2869 ! write (iout,*) "Calculating sidechain center iii",iii
2872 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2875 call sccenter(ires_old,iii,sccor)
2879 ! Start new residue.
2880 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2883 else if (ibeg.eq.1) then
2884 write (iout,*) "BEG ires",ires
2886 if (res.ne.'GLY' .and. res.ne. 'ACE') then
2889 nres_molec(molecule)=nres_molec(molecule)+1
2891 ires=ires-ishift+ishift1
2893 ! write (iout,*) "ishift",ishift," ires",ires,&
2894 ! " ires_old",ires_old
2896 else if (ibeg.eq.2) then
2898 ishift=-ires_old+ires-1 !!!!!
2899 ishift1=ishift1-1 !!!!!
2900 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2901 ires=ires-ishift+ishift1
2902 ! print *,ires,ishift,ishift1
2906 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2907 ires=ires-ishift+ishift1
2910 ! print *,'atom',ires,atom
2911 if (res.eq.'ACE' .or. res.eq.'NHE') then
2914 if (atom.eq.'CA '.or.atom.eq.'N ') then
2916 itype(ires,molecule)=rescode(ires,res,0,molecule)
2918 ! nres_molec(molecule)=nres_molec(molecule)+1
2921 itype(ires,molecule)=rescode(ires,res(2:3),0,molecule)
2922 ! nres_molec(molecule)=nres_molec(molecule)+1
2923 read (card(19:19),'(a1)') sugar
2924 isugar=sugarcode(sugar,ires)
2925 ! if (ibeg.eq.1) then
2929 ! print *,"ires=",ires,istype(ires)
2935 ires=ires-ishift+ishift1
2937 ! write (iout,*) "ires_old",ires_old," ires",ires
2938 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2941 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
2942 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2943 res.eq.'NHE'.and.atom(:2).eq.'HN') then
2944 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2945 ! print *,ires,ishift,ishift1
2946 ! write (iout,*) "backbone ",atom
2948 write (iout,'(2i3,2x,a,3f8.3)') &
2949 ires,itype(ires,1),res,(c(j,ires),j=1,3)
2952 nres_molec(molecule)=nres_molec(molecule)+1
2954 sccor(j,iii)=c(j,ires)
2956 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
2957 atom.eq."C2'" .or. atom.eq."C3'" &
2958 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
2959 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
2960 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
2961 ! print *,ires,ishift,ishift1
2965 ! sccor(j,iii)=c(j,ires)
2968 c(j,ires)=c(j,ires)+ccc(j)/5.0
2970 print *,counter,molecule
2971 if (counter.eq.5) then
2973 nres_molec(molecule)=nres_molec(molecule)+1
2976 ! sccor(j,iii)=c(j,ires)
2980 ! print *, "ATOM",atom(1:3)
2981 else if (atom.eq."C5'") then
2982 read (card(19:19),'(a1)') sugar
2983 isugar=sugarcode(sugar,ires)
2988 ! print *,ires,istype(ires)
2992 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
2993 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2994 nres_molec(molecule)=nres_molec(molecule)+1
2995 print *,"nres_molec(molecule)",nres_molec(molecule),ires
2999 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3001 else if ((atom.eq."C1'").and.unres_pdb) then
3003 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3004 ! write (*,*) card(23:27),ires,itype(ires,1)
3005 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
3006 atom.ne.'N' .and. atom.ne.'C' .and. &
3007 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
3008 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
3009 .and. atom.ne.'P '.and. &
3010 atom(1:1).ne.'H' .and. &
3011 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
3013 ! write (iout,*) "sidechain ",atom
3014 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3015 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3016 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3018 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3021 else if ((ions).and.(card(1:6).eq.'HETATM')) then
3022 if (firstion.eq.0) then
3026 dc(j,ires)=sccor(j,iii)
3029 call sccenter(ires,iii,sccor)
3032 read (card(12:16),*) atom
3033 print *,"HETATOM", atom
3034 read (card(18:20),'(a3)') res
3035 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3036 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3037 .or.(atom(1:2).eq.'K ')) &
3040 if (molecule.ne.5) molecprev=molecule
3042 nres_molec(molecule)=nres_molec(molecule)+1
3043 print *,"HERE",nres_molec(molecule)
3045 itype(ires,molecule)=rescode(ires,res,0,molecule)
3046 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3050 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3051 if (ires.eq.0) return
3052 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3055 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3056 .and.(.not.unres_pdb)) &
3057 nres_molec(molecule)=nres_molec(molecule)-2
3058 print *,'I have',nres, nres_molec(:)
3060 do k=1,4 ! ions are without dummy
3061 if (nres_molec(k).eq.0) cycle
3063 ! write (iout,*) i,itype(i,1)
3064 ! if (itype(i,1).eq.ntyp1) then
3065 ! write (iout,*) "dummy",i,itype(i,1)
3067 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3068 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3072 if (itype(i,k).eq.ntyp1_molec(k)) then
3073 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3074 if (itype(i+2,k).eq.0) then
3075 ! print *,"masz sieczke"
3077 if (itype(i+2,j).ne.0) then
3079 itype(i+1,j)=ntyp1_molec(j)
3080 nres_molec(k)=nres_molec(k)-1
3081 nres_molec(j)=nres_molec(j)+1
3087 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3088 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3089 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3090 ! if (unres_pdb) then
3091 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3092 ! print *,i,'tu dochodze'
3093 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3101 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3105 dcj=(c(j,i-2)-c(j,i-3))/2.0
3106 if (dcj.eq.0) dcj=1.23591524223
3111 else !itype(i+1,1).eq.ntyp1
3112 ! if (unres_pdb) then
3113 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3114 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3121 c(j,i)=c(j,i+1)-1.9d0*e2(j)
3125 dcj=(c(j,i+3)-c(j,i+2))/2.0
3126 if (dcj.eq.0) dcj=1.23591524223
3131 endif !itype(i+1,1).eq.ntyp1
3132 endif !itype.eq.ntyp1
3136 ! Calculate the CM of the last side chain.
3140 dc(j,ires)=sccor(j,iii)
3143 call sccenter(ires,iii,sccor)
3149 ! print *,"molecule",molecule
3150 if ((itype(nres,1).ne.10)) then
3152 if (molecule.eq.5) molecule=molecprev
3153 itype(nres,molecule)=ntyp1_molec(molecule)
3154 nres_molec(molecule)=nres_molec(molecule)+1
3155 ! if (unres_pdb) then
3156 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3157 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3164 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3168 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3169 c(j,nres)=c(j,nres-1)+dcj
3170 c(j,2*nres)=c(j,nres)
3174 ! print *,'I have',nres, nres_molec(:)
3176 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3178 if (nres.ne.nres0) then
3179 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3181 stop "Error nres value in WHAM input"
3184 !---------------------------------
3185 !el reallocate tables
3188 ! hfrag_alloc(j,i)=hfrag(j,i)
3191 ! bfrag_alloc(j,i)=bfrag(j,i)
3197 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3198 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3199 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3200 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3204 ! hfrag(j,i)=hfrag_alloc(j,i)
3209 ! bfrag(j,i)=bfrag_alloc(j,i)
3212 !el end reallocate tables
3213 !---------------------------------
3221 c(j,2*nres)=c(j,nres)
3224 if (itype(1,1).eq.ntyp1) then
3228 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3229 call refsys(2,3,4,e1,e2,e3,fail)
3236 c(j,1)=c(j,2)-1.9d0*e2(j)
3240 dcj=(c(j,4)-c(j,3))/2.0
3246 ! First lets assign correct dummy to correct type of chain
3248 if (itype(1,1).eq.ntyp1) then
3249 if (itype(2,1).eq.0) then
3251 if (itype(2,j).ne.0) then
3253 itype(1,j)=ntyp1_molec(j)
3254 nres_molec(1)=nres_molec(1)-1
3255 nres_molec(j)=nres_molec(j)+1
3262 print *,'I have',nres, nres_molec(:)
3264 ! Copy the coordinates to reference coordinates
3270 ! Calculate internal coordinates.
3272 write (iout,'(/a)') &
3273 "Cartesian coordinates of the reference structure"
3274 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3275 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3277 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3278 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3279 (c(j,ires+nres),j=1,3)
3282 ! znamy już nres wiec mozna alokowac tablice
3283 ! Calculate internal coordinates.
3284 if(me.eq.king.or..not.out1file)then
3285 write (iout,'(a)') &
3286 "Backbone and SC coordinates as read from the PDB"
3288 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
3289 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3290 (c(j,nres+ires),j=1,3)
3293 ! NOW LETS ROCK! SORTING
3294 allocate(c_temporary(3,2*nres))
3295 allocate(itype_temporary(nres,5))
3296 if (.not.allocated(molnum)) allocate(molnum(nres+1))
3297 if (.not.allocated(istype)) write(iout,*) &
3298 "SOMETHING WRONG WITH ISTYTPE"
3299 allocate(istype_temp(nres))
3300 itype_temporary(:,:)=0
3304 if (itype(i,k).ne.0) then
3306 c_temporary(j,seqalingbegin)=c(j,i)
3307 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3310 itype_temporary(seqalingbegin,k)=itype(i,k)
3311 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3312 istype_temp(seqalingbegin)=istype(i)
3313 molnum(seqalingbegin)=k
3314 seqalingbegin=seqalingbegin+1
3320 c(j,i)=c_temporary(j,i)
3325 itype(i,k)=itype_temporary(i,k)
3326 istype(i)=istype_temp(i)
3329 if (itype(1,1).eq.ntyp1) then
3333 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3334 call refsys(2,3,4,e1,e2,e3,fail)
3341 c(j,1)=c(j,2)-1.9d0*e2(j)
3345 dcj=(c(j,4)-c(j,3))/2.0
3353 write (iout,'(/a)') &
3354 "Cartesian coordinates of the reference structure after sorting"
3355 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
3356 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3358 write (iout,'(5(a3,1x),i3,3f8.3,5x,3f8.3)') &
3359 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3360 (c(j,ires+nres),j=1,3)
3364 ! print *,seqalingbegin,nres
3365 if(.not.allocated(vbld)) then
3366 allocate(vbld(2*nres))
3371 if(.not.allocated(vbld_inv)) then
3372 allocate(vbld_inv(2*nres))
3378 if(.not.allocated(theta)) then
3379 allocate(theta(nres+2))
3383 if(.not.allocated(phi)) allocate(phi(nres+2))
3384 if(.not.allocated(alph)) allocate(alph(nres+2))
3385 if(.not.allocated(omeg)) allocate(omeg(nres+2))
3386 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3387 if(.not.allocated(phiref)) allocate(phiref(nres+2))
3388 if(.not.allocated(costtab)) allocate(costtab(nres))
3389 if(.not.allocated(sinttab)) allocate(sinttab(nres))
3390 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3391 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3392 if(.not.allocated(xxref)) allocate(xxref(nres))
3393 if(.not.allocated(yyref)) allocate(yyref(nres))
3394 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3395 if(.not.allocated(dc_norm)) then
3396 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3397 allocate(dc_norm(3,0:2*nres+2))
3401 call int_from_cart(.true.,.false.)
3402 call sc_loc_geom(.false.)
3404 thetaref(i)=theta(i)
3414 dc(j,i)=c(j,i+1)-c(j,i)
3415 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3420 dc(j,i+nres)=c(j,i+nres)-c(j,i)
3421 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3423 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3427 ! Copy the coordinates to reference coordinates
3428 ! Splits to single chain if occurs
3432 ! cref(j,i,cou)=c(j,i)
3436 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3437 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
3438 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3439 !-----------------------------
3443 write (iout,*) "symetr", symetr
3446 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3448 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
3451 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3456 cref(j,i,cou)=c(j,i)
3457 cref(j,i+nres,cou)=c(j,i+nres)
3459 chain_rep(j,lll,kkk)=c(j,i)
3460 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
3464 write (iout,*) chain_length
3465 if (chain_length.eq.0) chain_length=nres
3467 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
3468 chain_rep(j,chain_length+nres,symetr) &
3469 =chain_rep(j,chain_length+nres,1)
3472 ! write (iout,*) "spraw lancuchy",chain_length,symetr
3474 ! do kkk=1,chain_length
3475 ! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
3479 ! makes copy of chains
3480 write (iout,*) "symetr", symetr
3485 if (symetr.gt.1) then
3492 write(iout,*) (tabperm(i,kkk),kkk=1,4)
3498 ! write (iout,*) i,icha
3499 do lll=1,chain_length
3501 if (cou.le.nres) then
3503 kupa=mod(lll,chain_length)
3504 iprzes=(kkk-1)*chain_length+lll
3505 if (kupa.eq.0) kupa=chain_length
3506 ! write (iout,*) "kupa", kupa
3507 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3508 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3515 !-koniec robienia kopii
3518 write (iout,*) "nowa struktura", nperm
3520 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3522 cref(3,i,kkk),cref(1,nres+i,kkk),&
3523 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3525 100 format (//' alpha-carbon coordinates ',&
3526 ' centroid coordinates'/ &
3527 ' ', 6X,'X',11X,'Y',11X,'Z', &
3528 10X,'X',11X,'Y',11X,'Z')
3529 110 format (a,'(',i3,')',6f12.5)
3535 bfrag(i,j)=bfrag(i,j)-ishift
3541 hfrag(i,j)=hfrag(i,j)-ishift
3547 end subroutine readpdb
3548 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3549 !-----------------------------------------------------------------------------
3551 !-----------------------------------------------------------------------------
3552 subroutine read_control
3566 use random, only: random_init
3567 ! implicit real*8 (a-h,o-z)
3568 ! include 'DIMENSIONS'
3570 use prng, only:prng_restart
3572 logical :: OKRandom!, prng_restart
3575 ! include 'COMMON.IOUNITS'
3576 ! include 'COMMON.TIME1'
3577 ! include 'COMMON.THREAD'
3578 ! include 'COMMON.SBRIDGE'
3579 ! include 'COMMON.CONTROL'
3580 ! include 'COMMON.MCM'
3581 ! include 'COMMON.MAP'
3582 ! include 'COMMON.HEADER'
3583 ! include 'COMMON.CSA'
3584 ! include 'COMMON.CHAIN'
3585 ! include 'COMMON.MUCA'
3586 ! include 'COMMON.MD'
3587 ! include 'COMMON.FFIELD'
3588 ! include 'COMMON.INTERACT'
3589 ! include 'COMMON.SETUP'
3590 !el integer :: KDIAG,ICORFL,IXDR
3591 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3592 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3593 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
3594 ! character(len=80) :: ucase
3595 character(len=640) :: controlcard
3597 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3603 read (INP,'(a)') titel
3604 call card_concat(controlcard,.true.)
3605 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3606 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3607 call reada(controlcard,'SEED',seed,0.0D0)
3608 call random_init(seed)
3609 ! Set up the time limit (caution! The time must be input in minutes!)
3610 read_cart=index(controlcard,'READ_CART').gt.0
3611 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3612 call readi(controlcard,'SYM',symetr,1)
3613 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3614 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3615 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3616 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3617 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3618 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3619 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3620 call reada(controlcard,'DRMS',drms,0.1D0)
3621 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3622 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3623 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3624 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3625 write (iout,'(a,f10.1)')'DRMS = ',drms
3626 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3627 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3629 call readi(controlcard,'NZ_START',nz_start,0)
3630 call readi(controlcard,'NZ_END',nz_end,0)
3631 call readi(controlcard,'IZ_SC',iz_sc,0)
3632 timlim=60.0D0*timlim
3633 safety = 60.0d0*safety
3636 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3637 !C SHIELD keyword sets if the shielding effect of side-chains is used
3638 !C 0 denots no shielding is used all peptide are equally despite the
3639 !C solvent accesible area
3640 !C 1 the newly introduced function
3641 !C 2 reseved for further possible developement
3642 call readi(controlcard,'SHIELD',shield_mode,0)
3643 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3644 write(iout,*) "shield_mode",shield_mode
3645 !C Varibles set size of box
3646 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3647 protein=index(controlcard,"PROTEIN").gt.0
3648 ions=index(controlcard,"IONS").gt.0
3649 nucleic=index(controlcard,"NUCLEIC").gt.0
3650 write (iout,*) "with_theta_constr ",with_theta_constr
3651 AFMlog=(index(controlcard,'AFM'))
3652 selfguide=(index(controlcard,'SELFGUIDE'))
3653 print *,'AFMlog',AFMlog,selfguide,"KUPA"
3654 call readi(controlcard,'GENCONSTR',genconstr,0)
3655 call reada(controlcard,'BOXX',boxxsize,100.0d0)
3656 call reada(controlcard,'BOXY',boxysize,100.0d0)
3657 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3658 call readi(controlcard,'TUBEMOD',tubemode,0)
3659 write (iout,*) TUBEmode,"TUBEMODE"
3660 if (TUBEmode.gt.0) then
3661 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3662 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3663 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3664 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3665 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3666 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3667 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3668 buftubebot=bordtubebot+tubebufthick
3669 buftubetop=bordtubetop-tubebufthick
3672 ! CUTOFFF ON ELECTROSTATICS
3673 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3674 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3675 write(iout,*) "R_CUT_ELE=",r_cut_ele
3676 ! Lipidic parameters
3677 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3678 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3679 if (lipthick.gt.0.0d0) then
3680 bordliptop=(boxzsize+lipthick)/2.0
3681 bordlipbot=bordliptop-lipthick
3682 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3683 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3684 buflipbot=bordlipbot+lipbufthick
3685 bufliptop=bordliptop-lipbufthick
3686 if ((lipbufthick*2.0d0).gt.lipthick) &
3687 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3688 endif !lipthick.gt.0
3689 write(iout,*) "bordliptop=",bordliptop
3690 write(iout,*) "bordlipbot=",bordlipbot
3691 write(iout,*) "bufliptop=",bufliptop
3692 write(iout,*) "buflipbot=",buflipbot
3693 write (iout,*) "SHIELD MODE",shield_mode
3695 !C-------------------------
3696 minim=(index(controlcard,'MINIMIZE').gt.0)
3697 dccart=(index(controlcard,'CART').gt.0)
3698 overlapsc=(index(controlcard,'OVERLAP').gt.0)
3699 overlapsc=.not.overlapsc
3700 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3701 searchsc=.not.searchsc
3702 sideadd=(index(controlcard,'SIDEADD').gt.0)
3703 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3704 outpdb=(index(controlcard,'PDBOUT').gt.0)
3705 outmol2=(index(controlcard,'MOL2OUT').gt.0)
3706 pdbref=(index(controlcard,'PDBREF').gt.0)
3707 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3708 indpdb=index(controlcard,'PDBSTART')
3709 extconf=(index(controlcard,'EXTCONF').gt.0)
3710 call readi(controlcard,'IPRINT',iprint,0)
3711 call readi(controlcard,'MAXGEN',maxgen,10000)
3712 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3713 call readi(controlcard,"KDIAG",kdiag,0)
3714 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3715 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3716 write (iout,*) "RESCALE_MODE",rescale_mode
3717 split_ene=index(controlcard,'SPLIT_ENE').gt.0
3718 if (index(controlcard,'REGULAR').gt.0.0D0) then
3719 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3723 if (index(controlcard,'CHECKGRAD').gt.0) then
3725 if (index(controlcard,'CART').gt.0) then
3727 elseif (index(controlcard,'CARINT').gt.0) then
3732 elseif (index(controlcard,'THREAD').gt.0) then
3734 call readi(controlcard,'THREAD',nthread,0)
3735 if (nthread.gt.0) then
3736 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3739 write (iout,'(a)')'A number has to follow the THREAD keyword.'
3740 stop 'Error termination in Read_Control.'
3742 else if (index(controlcard,'MCMA').gt.0) then
3744 else if (index(controlcard,'MCEE').gt.0) then
3746 else if (index(controlcard,'MULTCONF').gt.0) then
3748 else if (index(controlcard,'MAP').gt.0) then
3750 call readi(controlcard,'MAP',nmap,0)
3751 else if (index(controlcard,'CSA').gt.0) then
3753 !rc else if (index(controlcard,'ZSCORE').gt.0) then
3755 !rc ZSCORE is rm from UNRES, modecalc=9 is available
3758 !fcm else if (index(controlcard,'MCMF').gt.0) then
3760 else if (index(controlcard,'SOFTREG').gt.0) then
3762 else if (index(controlcard,'CHECK_BOND').gt.0) then
3764 else if (index(controlcard,'TEST').gt.0) then
3766 else if (index(controlcard,'MD').gt.0) then
3768 else if (index(controlcard,'RE ').gt.0) then
3772 lmuca=index(controlcard,'MUCA').gt.0
3773 call readi(controlcard,'MUCADYN',mucadyn,0)
3774 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3775 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3777 write (iout,*) 'MUCADYN=',mucadyn
3778 write (iout,*) 'MUCASMOOTH=',muca_smooth
3781 iscode=index(controlcard,'ONE_LETTER')
3782 indphi=index(controlcard,'PHI')
3783 indback=index(controlcard,'BACK')
3784 iranconf=index(controlcard,'RAND_CONF')
3785 i2ndstr=index(controlcard,'USE_SEC_PRED')
3786 gradout=index(controlcard,'GRADOUT').gt.0
3787 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3788 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3789 if (me.eq.king .or. .not.out1file ) &
3790 write (iout,*) "DISTCHAINMAX",distchainmax
3792 if(me.eq.king.or..not.out1file) &
3793 write (iout,'(2a)') diagmeth(kdiag),&
3794 ' routine used to diagonalize matrices.'
3795 if (shield_mode.gt.0) then
3797 !C VSolvSphere the volume of solving sphere
3799 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3800 !C there will be no distinction between proline peptide group and normal peptide
3801 !C group in case of shielding parameters
3802 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3803 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3804 write (iout,*) VSolvSphere,VSolvSphere_div
3805 !C long axis of side chain
3807 long_r_sidechain(i)=vbldsc0(1,i)
3808 short_r_sidechain(i)=sigma0(i)
3809 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3815 end subroutine read_control
3816 !-----------------------------------------------------------------------------
3817 subroutine read_REMDpar
3819 ! Read REMD settings
3826 use control_data, only:out1file
3827 ! implicit real*8 (a-h,o-z)
3828 ! include 'DIMENSIONS'
3829 ! include 'COMMON.IOUNITS'
3830 ! include 'COMMON.TIME1'
3831 ! include 'COMMON.MD'
3834 !el include 'COMMON.LANGEVIN'
3836 !el include 'COMMON.LANGEVIN.lang0'
3838 ! include 'COMMON.INTERACT'
3839 ! include 'COMMON.NAMES'
3840 ! include 'COMMON.GEO'
3841 ! include 'COMMON.REMD'
3842 ! include 'COMMON.CONTROL'
3843 ! include 'COMMON.SETUP'
3844 ! character(len=80) :: ucase
3845 character(len=320) :: controlcard
3846 character(len=3200) :: controlcard1
3847 integer :: iremd_m_total
3850 ! real(kind=8) :: var,ene
3852 if(me.eq.king.or..not.out1file) &
3853 write (iout,*) "REMD setup"
3855 call card_concat(controlcard,.true.)
3856 call readi(controlcard,"NREP",nrep,3)
3857 call readi(controlcard,"NSTEX",nstex,1000)
3858 call reada(controlcard,"RETMIN",retmin,10.0d0)
3859 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3860 mremdsync=(index(controlcard,'SYNC').gt.0)
3861 call readi(controlcard,"NSYN",i_sync_step,100)
3862 restart1file=(index(controlcard,'REST1FILE').gt.0)
3863 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3864 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3865 if(max_cache_traj_use.gt.max_cache_traj) &
3866 max_cache_traj_use=max_cache_traj
3867 if(me.eq.king.or..not.out1file) then
3868 !d if (traj1file) then
3869 !rc caching is in testing - NTWX is not ignored
3870 !d write (iout,*) "NTWX value is ignored"
3871 !d write (iout,*) " trajectory is stored to one file by master"
3872 !d write (iout,*) " before exchange at NSTEX intervals"
3874 write (iout,*) "NREP= ",nrep
3875 write (iout,*) "NSTEX= ",nstex
3876 write (iout,*) "SYNC= ",mremdsync
3877 write (iout,*) "NSYN= ",i_sync_step
3878 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3881 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3882 if (index(controlcard,'TLIST').gt.0) then
3884 call card_concat(controlcard1,.true.)
3885 read(controlcard1,*) (remd_t(i),i=1,nrep)
3886 if(me.eq.king.or..not.out1file) &
3887 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3890 if (index(controlcard,'MLIST').gt.0) then
3892 call card_concat(controlcard1,.true.)
3893 read(controlcard1,*) (remd_m(i),i=1,nrep)
3894 if(me.eq.king.or..not.out1file) then
3895 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3898 iremd_m_total=iremd_m_total+remd_m(i)
3900 write (iout,*) 'Total number of replicas ',iremd_m_total
3903 if(me.eq.king.or..not.out1file) &
3904 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3906 end subroutine read_REMDpar
3907 !-----------------------------------------------------------------------------
3908 subroutine read_MDpar
3912 use control_data, only: r_cut,rlamb,out1file
3914 use geometry_data, only: pi
3916 ! implicit real*8 (a-h,o-z)
3917 ! include 'DIMENSIONS'
3918 ! include 'COMMON.IOUNITS'
3919 ! include 'COMMON.TIME1'
3920 ! include 'COMMON.MD'
3923 !el include 'COMMON.LANGEVIN'
3925 !el include 'COMMON.LANGEVIN.lang0'
3927 ! include 'COMMON.INTERACT'
3928 ! include 'COMMON.NAMES'
3929 ! include 'COMMON.GEO'
3930 ! include 'COMMON.SETUP'
3931 ! include 'COMMON.CONTROL'
3932 ! include 'COMMON.SPLITELE'
3933 ! character(len=80) :: ucase
3934 character(len=320) :: controlcard
3939 call card_concat(controlcard,.true.)
3940 call readi(controlcard,"NSTEP",n_timestep,1000000)
3941 call readi(controlcard,"NTWE",ntwe,100)
3942 call readi(controlcard,"NTWX",ntwx,1000)
3943 call reada(controlcard,"DT",d_time,1.0d-1)
3944 call reada(controlcard,"DVMAX",dvmax,2.0d1)
3945 call reada(controlcard,"DAMAX",damax,1.0d1)
3946 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3947 call readi(controlcard,"LANG",lang,0)
3948 RESPA = index(controlcard,"RESPA") .gt. 0
3949 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3950 ntime_split0=ntime_split
3951 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3952 ntime_split0=ntime_split
3953 call reada(controlcard,"R_CUT",r_cut,2.0d0)
3954 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3955 rest = index(controlcard,"REST").gt.0
3956 tbf = index(controlcard,"TBF").gt.0
3957 usampl = index(controlcard,"USAMPL").gt.0
3958 mdpdb = index(controlcard,"MDPDB").gt.0
3959 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3960 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3961 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3962 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3963 if (count_reset_moment.eq.0) count_reset_moment=1000000000
3964 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3965 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3966 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3967 if (count_reset_vel.eq.0) count_reset_vel=1000000000
3968 large = index(controlcard,"LARGE").gt.0
3969 print_compon = index(controlcard,"PRINT_COMPON").gt.0
3970 rattle = index(controlcard,"RATTLE").gt.0
3971 ! if performing umbrella sampling, fragments constrained are read from the fragment file
3977 if(me.eq.king.or..not.out1file) then
3979 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3981 write (iout,'(a)') "The units are:"
3982 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3983 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3984 " acceleration: angstrom/(48.9 fs)**2"
3985 write (iout,'(a)') "energy: kcal/mol, temperature: K"
3987 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3988 write (iout,'(a60,f10.5,a)') &
3989 "Initial time step of numerical integration:",d_time,&
3991 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3993 write (iout,'(2a,i4,a)') &
3994 "A-MTS algorithm used; initial time step for fast-varying",&
3995 " short-range forces split into",ntime_split," steps."
3996 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3997 r_cut," lambda",rlamb
3999 write (iout,'(2a,f10.5)') &
4000 "Maximum acceleration threshold to reduce the time step",&
4001 "/increase split number:",damax
4002 write (iout,'(2a,f10.5)') &
4003 "Maximum predicted energy drift to reduce the timestep",&
4004 "/increase split number:",edriftmax
4005 write (iout,'(a60,f10.5)') &
4006 "Maximum velocity threshold to reduce velocities:",dvmax
4007 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
4008 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
4009 if (rattle) write (iout,'(a60)') &
4010 "Rattle algorithm used to constrain the virtual bonds"
4014 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
4015 call reada(controlcard,"RWAT",rwat,1.4d0)
4016 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
4017 surfarea=index(controlcard,"SURFAREA").gt.0
4018 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4019 if(me.eq.king.or..not.out1file)then
4020 write (iout,'(/a,$)') "Langevin dynamics calculation"
4022 write (iout,'(a/)') &
4023 " with direct integration of Langevin equations"
4024 else if (lang.eq.2) then
4025 write (iout,'(a/)') " with TINKER stochasic MD integrator"
4026 else if (lang.eq.3) then
4027 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4028 else if (lang.eq.4) then
4029 write (iout,'(a/)') " in overdamped mode"
4031 write (iout,'(//a,i5)') &
4032 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4035 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4036 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4037 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4038 write (iout,'(a60,f10.5)') &
4039 "Scaling factor of the friction forces:",scal_fric
4040 if (surfarea) write (iout,'(2a,i10,a)') &
4041 "Friction coefficients will be scaled by solvent-accessible",&
4042 " surface area every",reset_fricmat," steps."
4044 ! Calculate friction coefficients and bounds of stochastic forces
4045 eta=6*pi*cPoise*etawat
4046 if(me.eq.king.or..not.out1file) &
4047 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4050 do j=1,5 !types of molecules
4051 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4052 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4054 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4055 do j=1,5 !types of molecules
4057 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4058 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4062 if(me.eq.king.or..not.out1file)then
4063 write (iout,'(/2a/)') &
4064 "Radii of site types and friction coefficients and std's of",&
4065 " stochastic forces of fully exposed sites"
4066 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4068 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4069 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4073 if(me.eq.king.or..not.out1file)then
4074 write (iout,'(a)') "Berendsen bath calculation"
4075 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4076 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4078 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4079 count_reset_moment," steps"
4081 write (iout,'(a,i10,a)') &
4082 "Velocities will be reset at random every",count_reset_vel,&
4086 if(me.eq.king.or..not.out1file) &
4087 write (iout,'(a31)') "Microcanonical mode calculation"
4089 if(me.eq.king.or..not.out1file)then
4090 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4092 write(iout,*) "MD running with constraints."
4093 write(iout,*) "Equilibration time ", eq_time, " mtus."
4094 write(iout,*) "Constraining ", nfrag," fragments."
4095 write(iout,*) "Length of each fragment, weight and q0:"
4097 write (iout,*) "Set of restraints #",iset
4099 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4100 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4102 write(iout,*) "constraints between ", npair, "fragments."
4103 write(iout,*) "constraint pairs, weights and q0:"
4105 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4106 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4108 write(iout,*) "angle constraints within ", nfrag_back,&
4109 "backbone fragments."
4110 write(iout,*) "fragment, weights:"
4112 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4113 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4114 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4117 iset=mod(kolor,nset)+1
4120 if(me.eq.king.or..not.out1file) &
4121 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4123 end subroutine read_MDpar
4124 !-----------------------------------------------------------------------------
4128 ! implicit real*8 (a-h,o-z)
4129 ! include 'DIMENSIONS'
4130 ! include 'COMMON.MAP'
4131 ! include 'COMMON.IOUNITS'
4132 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4133 character(len=80) :: mapcard !,ucase
4136 ! real(kind=8) :: var,ene
4139 read (inp,'(a)') mapcard
4140 mapcard=ucase(mapcard)
4141 if (index(mapcard,'PHI').gt.0) then
4143 else if (index(mapcard,'THE').gt.0) then
4145 else if (index(mapcard,'ALP').gt.0) then
4147 else if (index(mapcard,'OME').gt.0) then
4150 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4151 stop 'Error - illegal variable spec in MAP card.'
4153 call readi (mapcard,'RES1',res1(imap),0)
4154 call readi (mapcard,'RES2',res2(imap),0)
4155 if (res1(imap).eq.0) then
4156 res1(imap)=res2(imap)
4157 else if (res2(imap).eq.0) then
4158 res2(imap)=res1(imap)
4160 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4161 write (iout,'(a)') &
4162 'Error - illegal definition of variable group in MAP.'
4163 stop 'Error - illegal definition of variable group in MAP.'
4165 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4166 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4167 call readi(mapcard,'NSTEP',nstep(imap),0)
4168 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4169 write (iout,'(a)') &
4170 'Illegal boundary and/or step size specification in MAP.'
4171 stop 'Illegal boundary and/or step size specification in MAP.'
4175 end subroutine map_read
4176 !-----------------------------------------------------------------------------
4179 use control_data, only: vdisulf
4181 ! implicit real*8 (a-h,o-z)
4182 ! include 'DIMENSIONS'
4183 ! include 'COMMON.IOUNITS'
4184 ! include 'COMMON.GEO'
4185 ! include 'COMMON.CSA'
4186 ! include 'COMMON.BANK'
4187 ! include 'COMMON.CONTROL'
4188 ! character(len=80) :: ucase
4189 character(len=620) :: mcmcard
4191 ! integer :: ntf,ik,iw_pdb
4192 ! real(kind=8) :: var,ene
4194 call card_concat(mcmcard,.true.)
4196 call readi(mcmcard,'NCONF',nconf,50)
4197 call readi(mcmcard,'NADD',nadd,0)
4198 call readi(mcmcard,'JSTART',jstart,1)
4199 call readi(mcmcard,'JEND',jend,1)
4200 call readi(mcmcard,'NSTMAX',nstmax,500000)
4201 call readi(mcmcard,'N0',n0,1)
4202 call readi(mcmcard,'N1',n1,6)
4203 call readi(mcmcard,'N2',n2,4)
4204 call readi(mcmcard,'N3',n3,0)
4205 call readi(mcmcard,'N4',n4,0)
4206 call readi(mcmcard,'N5',n5,0)
4207 call readi(mcmcard,'N6',n6,10)
4208 call readi(mcmcard,'N7',n7,0)
4209 call readi(mcmcard,'N8',n8,0)
4210 call readi(mcmcard,'N9',n9,0)
4211 call readi(mcmcard,'N14',n14,0)
4212 call readi(mcmcard,'N15',n15,0)
4213 call readi(mcmcard,'N16',n16,0)
4214 call readi(mcmcard,'N17',n17,0)
4215 call readi(mcmcard,'N18',n18,0)
4217 vdisulf=(index(mcmcard,'DYNSS').gt.0)
4219 call readi(mcmcard,'NDIFF',ndiff,2)
4220 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4221 call readi(mcmcard,'IS1',is1,1)
4222 call readi(mcmcard,'IS2',is2,8)
4223 call readi(mcmcard,'NRAN0',nran0,4)
4224 call readi(mcmcard,'NRAN1',nran1,2)
4225 call readi(mcmcard,'IRR',irr,1)
4226 call readi(mcmcard,'NSEED',nseed,20)
4227 call readi(mcmcard,'NTOTAL',ntotal,10000)
4228 call reada(mcmcard,'CUT1',cut1,2.0d0)
4229 call reada(mcmcard,'CUT2',cut2,5.0d0)
4230 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4231 call readi(mcmcard,'ICMAX',icmax,3)
4232 call readi(mcmcard,'IRESTART',irestart,0)
4233 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
4236 call reada(mcmcard,'DELE',dele,20.0d0)
4237 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4238 call readi(mcmcard,'IREF',iref,0)
4239 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4240 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4241 call readi(mcmcard,'NCONF_IN',nconf_in,0)
4242 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4243 write (iout,*) "NCONF_IN",nconf_in
4245 end subroutine csaread
4246 !-----------------------------------------------------------------------------
4250 use control_data, only: MaxMoveType
4253 ! implicit real*8 (a-h,o-z)
4254 ! include 'DIMENSIONS'
4255 ! include 'COMMON.MCM'
4256 ! include 'COMMON.MCE'
4257 ! include 'COMMON.IOUNITS'
4258 ! character(len=80) :: ucase
4259 character(len=320) :: mcmcard
4262 ! real(kind=8) :: var,ene
4264 call card_concat(mcmcard,.true.)
4265 call readi(mcmcard,'MAXACC',maxacc,100)
4266 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4267 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4268 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4269 call readi(mcmcard,'MAXREPM',maxrepm,200)
4270 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4271 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4272 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4273 call reada(mcmcard,'E_UP',e_up,5.0D0)
4274 call reada(mcmcard,'DELTE',delte,0.1D0)
4275 call readi(mcmcard,'NSWEEP',nsweep,5)
4276 call readi(mcmcard,'NSTEPH',nsteph,0)
4277 call readi(mcmcard,'NSTEPC',nstepc,0)
4278 call reada(mcmcard,'TMIN',tmin,298.0D0)
4279 call reada(mcmcard,'TMAX',tmax,298.0D0)
4280 call readi(mcmcard,'NWINDOW',nwindow,0)
4281 call readi(mcmcard,'PRINT_MC',print_mc,0)
4282 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4283 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4284 ent_read=(index(mcmcard,'ENT_READ').gt.0)
4285 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4286 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4287 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4288 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4289 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4290 if (nwindow.gt.0) then
4291 allocate(winstart(nwindow)) !!el (maxres)
4292 allocate(winend(nwindow)) !!el
4293 allocate(winlen(nwindow)) !!el
4294 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4296 winlen(i)=winend(i)-winstart(i)+1
4299 if (tmax.lt.tmin) tmax=tmin
4300 if (tmax.eq.tmin) then
4304 if (nstepc.gt.0 .and. nsteph.gt.0) then
4305 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4306 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4308 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4309 ! Probabilities of different move types
4310 sumpro_type(0)=0.0D0
4311 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4312 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4313 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4314 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
4315 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4316 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4317 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4319 print *,'i',i,' sumprotype',sumpro_type(i)
4320 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4321 print *,'i',i,' sumprotype',sumpro_type(i)
4324 end subroutine mcmread
4325 !-----------------------------------------------------------------------------
4326 subroutine read_minim
4329 ! implicit real*8 (a-h,o-z)
4330 ! include 'DIMENSIONS'
4331 ! include 'COMMON.MINIM'
4332 ! include 'COMMON.IOUNITS'
4333 ! character(len=80) :: ucase
4334 character(len=320) :: minimcard
4336 ! integer :: ntf,ik,iw_pdb
4337 ! real(kind=8) :: var,ene
4339 call card_concat(minimcard,.true.)
4340 call readi(minimcard,'MAXMIN',maxmin,2000)
4341 call readi(minimcard,'MAXFUN',maxfun,5000)
4342 call readi(minimcard,'MINMIN',minmin,maxmin)
4343 call readi(minimcard,'MINFUN',minfun,maxmin)
4344 call reada(minimcard,'TOLF',tolf,1.0D-2)
4345 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4346 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4347 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4348 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4349 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4350 'Options in energy minimization:'
4351 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4352 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4353 'MinMin:',MinMin,' MinFun:',MinFun,&
4354 ' TolF:',TolF,' RTolF:',RTolF
4356 end subroutine read_minim
4357 !-----------------------------------------------------------------------------
4358 subroutine openunits
4360 use MD_data, only: usampl
4363 use control_data, only:out1file
4364 use control, only: getenv_loc
4365 ! implicit real*8 (a-h,o-z)
4366 ! include 'DIMENSIONS'
4369 character(len=16) :: form,nodename
4370 integer :: nodelen,ierror,npos
4372 ! include 'COMMON.SETUP'
4373 ! include 'COMMON.IOUNITS'
4374 ! include 'COMMON.MD'
4375 ! include 'COMMON.CONTROL'
4376 integer :: lenpre,lenpot,lentmp !,ilen
4378 character(len=3) :: out1file_text !,ucase
4379 character(len=3) :: ll
4382 ! integer :: ntf,ik,iw_pdb
4383 ! real(kind=8) :: var,ene
4385 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4386 call getenv_loc("PREFIX",prefix)
4388 call getenv_loc("POT",pot)
4389 call getenv_loc("DIRTMP",tmpdir)
4390 call getenv_loc("CURDIR",curdir)
4391 call getenv_loc("OUT1FILE",out1file_text)
4392 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4393 out1file_text=ucase(out1file_text)
4394 if (out1file_text(1:1).eq."Y") then
4397 out1file=fg_rank.gt.0
4402 if (lentmp.gt.0) then
4403 write (*,'(80(1h!))')
4404 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
4405 write (*,'(80(1h!))')
4406 write (*,*)"All output files will be on node /tmp directory."
4408 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4409 if (me.eq.king) then
4410 write (*,*) "The master node is ",nodename
4411 else if (fg_rank.eq.0) then
4412 write (*,*) "I am the CG slave node ",nodename
4414 write (*,*) "I am the FG slave node ",nodename
4417 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4418 lenpre = lentmp+lenpre+1
4420 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4421 ! Get the names and open the input files
4422 #if defined(WINIFL) || defined(WINPGI)
4423 open(1,file=pref_orig(:ilen(pref_orig))// &
4424 '.inp',status='old',readonly,shared)
4425 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4426 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4427 ! Get parameter filenames and open the parameter files.
4428 call getenv_loc('BONDPAR',bondname)
4429 open (ibond,file=bondname,status='old',readonly,shared)
4430 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4431 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
4432 call getenv_loc('THETPAR',thetname)
4433 open (ithep,file=thetname,status='old',readonly,shared)
4434 call getenv_loc('ROTPAR',rotname)
4435 open (irotam,file=rotname,status='old',readonly,shared)
4436 call getenv_loc('TORPAR',torname)
4437 open (itorp,file=torname,status='old',readonly,shared)
4438 call getenv_loc('TORDPAR',tordname)
4439 open (itordp,file=tordname,status='old',readonly,shared)
4440 call getenv_loc('FOURIER',fouriername)
4441 open (ifourier,file=fouriername,status='old',readonly,shared)
4442 call getenv_loc('ELEPAR',elename)
4443 open (ielep,file=elename,status='old',readonly,shared)
4444 call getenv_loc('SIDEPAR',sidename)
4445 open (isidep,file=sidename,status='old',readonly,shared)
4447 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4448 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
4449 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4450 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
4451 call getenv_loc('TORPAR_NUCL',torname_nucl)
4452 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
4453 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4454 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
4455 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4456 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
4459 #elif (defined CRAY) || (defined AIX)
4460 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4462 ! print *,"Processor",myrank," opened file 1"
4463 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4464 ! print *,"Processor",myrank," opened file 9"
4465 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4466 ! Get parameter filenames and open the parameter files.
4467 call getenv_loc('BONDPAR',bondname)
4468 open (ibond,file=bondname,status='old',action='read')
4469 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4470 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4472 ! print *,"Processor",myrank," opened file IBOND"
4473 call getenv_loc('THETPAR',thetname)
4474 open (ithep,file=thetname,status='old',action='read')
4475 ! print *,"Processor",myrank," opened file ITHEP"
4476 call getenv_loc('ROTPAR',rotname)
4477 open (irotam,file=rotname,status='old',action='read')
4478 ! print *,"Processor",myrank," opened file IROTAM"
4479 call getenv_loc('TORPAR',torname)
4480 open (itorp,file=torname,status='old',action='read')
4481 ! print *,"Processor",myrank," opened file ITORP"
4482 call getenv_loc('TORDPAR',tordname)
4483 open (itordp,file=tordname,status='old',action='read')
4484 ! print *,"Processor",myrank," opened file ITORDP"
4485 call getenv_loc('SCCORPAR',sccorname)
4486 open (isccor,file=sccorname,status='old',action='read')
4487 ! print *,"Processor",myrank," opened file ISCCOR"
4488 call getenv_loc('FOURIER',fouriername)
4489 open (ifourier,file=fouriername,status='old',action='read')
4490 ! print *,"Processor",myrank," opened file IFOURIER"
4491 call getenv_loc('ELEPAR',elename)
4492 open (ielep,file=elename,status='old',action='read')
4493 ! print *,"Processor",myrank," opened file IELEP"
4494 call getenv_loc('SIDEPAR',sidename)
4495 open (isidep,file=sidename,status='old',action='read')
4497 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4498 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4499 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4500 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4501 call getenv_loc('TORPAR_NUCL',torname_nucl)
4502 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4503 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4504 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4505 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4506 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4508 call getenv_loc('LIPTRANPAR',liptranname)
4509 open (iliptranpar,file=liptranname,status='old',action='read')
4510 call getenv_loc('TUBEPAR',tubename)
4511 open (itube,file=tubename,status='old',action='read')
4512 call getenv_loc('IONPAR',ionname)
4513 open (iion,file=ionname,status='old',action='read')
4515 ! print *,"Processor",myrank," opened file ISIDEP"
4516 ! print *,"Processor",myrank," opened parameter files"
4518 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
4519 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4520 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4521 ! Get parameter filenames and open the parameter files.
4522 call getenv_loc('BONDPAR',bondname)
4523 open (ibond,file=bondname,status='old')
4524 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4525 open (ibond_nucl,file=bondname_nucl,status='old')
4527 call getenv_loc('THETPAR',thetname)
4528 open (ithep,file=thetname,status='old')
4529 call getenv_loc('ROTPAR',rotname)
4530 open (irotam,file=rotname,status='old')
4531 call getenv_loc('TORPAR',torname)
4532 open (itorp,file=torname,status='old')
4533 call getenv_loc('TORDPAR',tordname)
4534 open (itordp,file=tordname,status='old')
4535 call getenv_loc('SCCORPAR',sccorname)
4536 open (isccor,file=sccorname,status='old')
4537 call getenv_loc('FOURIER',fouriername)
4538 open (ifourier,file=fouriername,status='old')
4539 call getenv_loc('ELEPAR',elename)
4540 open (ielep,file=elename,status='old')
4541 call getenv_loc('SIDEPAR',sidename)
4542 open (isidep,file=sidename,status='old')
4544 open (ithep_nucl,file=thetname_nucl,status='old')
4545 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4546 open (irotam_nucl,file=rotname_nucl,status='old')
4547 call getenv_loc('TORPAR_NUCL',torname_nucl)
4548 open (itorp_nucl,file=torname_nucl,status='old')
4549 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4550 open (itordp_nucl,file=tordname_nucl,status='old')
4551 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4552 open (isidep_nucl,file=sidename_nucl,status='old')
4554 call getenv_loc('LIPTRANPAR',liptranname)
4555 open (iliptranpar,file=liptranname,status='old')
4556 call getenv_loc('TUBEPAR',tubename)
4557 open (itube,file=tubename,status='old')
4558 call getenv_loc('IONPAR',ionname)
4559 open (iion,file=ionname,status='old')
4561 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4563 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4564 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4565 ! Get parameter filenames and open the parameter files.
4566 call getenv_loc('BONDPAR',bondname)
4567 open (ibond,file=bondname,status='old',action='read')
4568 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4569 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4570 call getenv_loc('THETPAR',thetname)
4571 open (ithep,file=thetname,status='old',action='read')
4572 call getenv_loc('ROTPAR',rotname)
4573 open (irotam,file=rotname,status='old',action='read')
4574 call getenv_loc('TORPAR',torname)
4575 open (itorp,file=torname,status='old',action='read')
4576 call getenv_loc('TORDPAR',tordname)
4577 open (itordp,file=tordname,status='old',action='read')
4578 call getenv_loc('SCCORPAR',sccorname)
4579 open (isccor,file=sccorname,status='old',action='read')
4581 call getenv_loc('THETPARPDB',thetname_pdb)
4582 print *,"thetname_pdb ",thetname_pdb
4583 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4584 print *,ithep_pdb," opened"
4586 call getenv_loc('FOURIER',fouriername)
4587 open (ifourier,file=fouriername,status='old',readonly)
4588 call getenv_loc('ELEPAR',elename)
4589 open (ielep,file=elename,status='old',readonly)
4590 call getenv_loc('SIDEPAR',sidename)
4591 open (isidep,file=sidename,status='old',readonly)
4593 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4594 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4595 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4596 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4597 call getenv_loc('TORPAR_NUCL',torname_nucl)
4598 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4599 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4600 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4601 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4602 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4603 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
4604 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
4605 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
4606 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
4607 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
4608 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
4609 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
4610 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
4613 call getenv_loc('LIPTRANPAR',liptranname)
4614 open (iliptranpar,file=liptranname,status='old',action='read')
4615 call getenv_loc('TUBEPAR',tubename)
4616 open (itube,file=tubename,status='old',action='read')
4617 call getenv_loc('IONPAR',ionname)
4618 open (iion,file=ionname,status='old',action='read')
4621 call getenv_loc('ROTPARPDB',rotname_pdb)
4622 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4625 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
4626 #if defined(WINIFL) || defined(WINPGI)
4627 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
4628 #elif (defined CRAY) || (defined AIX)
4629 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4631 open (iscpp_nucl,file=scpname_nucl,status='old')
4633 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4638 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4639 ! Use -DOLDSCP to use hard-coded constants instead.
4641 call getenv_loc('SCPPAR',scpname)
4642 #if defined(WINIFL) || defined(WINPGI)
4643 open (iscpp,file=scpname,status='old',readonly,shared)
4644 #elif (defined CRAY) || (defined AIX)
4645 open (iscpp,file=scpname,status='old',action='read')
4647 open (iscpp,file=scpname,status='old')
4649 open (iscpp,file=scpname,status='old',action='read')
4652 call getenv_loc('PATTERN',patname)
4653 #if defined(WINIFL) || defined(WINPGI)
4654 open (icbase,file=patname,status='old',readonly,shared)
4655 #elif (defined CRAY) || (defined AIX)
4656 open (icbase,file=patname,status='old',action='read')
4658 open (icbase,file=patname,status='old')
4660 open (icbase,file=patname,status='old',action='read')
4663 ! Open output file only for CG processes
4664 ! print *,"Processor",myrank," fg_rank",fg_rank
4665 if (fg_rank.eq.0) then
4667 if (nodes.eq.1) then
4670 npos = dlog10(dfloat(nodes-1))+1
4672 if (npos.lt.3) npos=3
4673 write (liczba,'(i1)') npos
4674 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4676 write (liczba,form) me
4677 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4678 liczba(:ilen(liczba))
4679 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4681 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4683 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4684 liczba(:ilen(liczba))//'.mol2'
4685 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4686 liczba(:ilen(liczba))//'.stat'
4688 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4689 //liczba(:ilen(liczba))//'.stat')
4690 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4693 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4694 liczba(:ilen(liczba))//'.const'
4699 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4700 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4701 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4702 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4703 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4705 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4707 rest2name=prefix(:ilen(prefix))//'.rst'
4709 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4712 #if defined(AIX) || defined(PGI)
4713 if (me.eq.king .or. .not. out1file) &
4714 open(iout,file=outname,status='unknown')
4716 if (fg_rank.gt.0) then
4717 write (liczba,'(i3.3)') myrank/nfgtasks
4718 write (ll,'(bz,i3.3)') fg_rank
4719 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4724 open(igeom,file=intname,status='unknown',position='append')
4725 open(ipdb,file=pdbname,status='unknown')
4726 open(imol2,file=mol2name,status='unknown')
4727 open(istat,file=statname,status='unknown',position='append')
4729 !1out open(iout,file=outname,status='unknown')
4732 if (me.eq.king .or. .not.out1file) &
4733 open(iout,file=outname,status='unknown')
4735 if (fg_rank.gt.0) then
4736 write (liczba,'(i3.3)') myrank/nfgtasks
4737 write (ll,'(bz,i3.3)') fg_rank
4738 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4743 open(igeom,file=intname,status='unknown',access='append')
4744 open(ipdb,file=pdbname,status='unknown')
4745 open(imol2,file=mol2name,status='unknown')
4746 open(istat,file=statname,status='unknown',access='append')
4748 !1out open(iout,file=outname,status='unknown')
4751 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4752 csa_seed=prefix(:lenpre)//'.CSA.seed'
4753 csa_history=prefix(:lenpre)//'.CSA.history'
4754 csa_bank=prefix(:lenpre)//'.CSA.bank'
4755 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4756 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4757 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4758 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4759 csa_int=prefix(:lenpre)//'.int'
4760 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4761 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4762 csa_in=prefix(:lenpre)//'.CSA.in'
4763 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4766 write (iout,'(80(1h-))')
4767 write (iout,'(30x,a)') "FILE ASSIGNMENT"
4768 write (iout,'(80(1h-))')
4769 write (iout,*) "Input file : ",&
4770 pref_orig(:ilen(pref_orig))//'.inp'
4771 write (iout,*) "Output file : ",&
4772 outname(:ilen(outname))
4774 write (iout,*) "Sidechain potential file : ",&
4775 sidename(:ilen(sidename))
4777 write (iout,*) "SCp potential file : ",&
4778 scpname(:ilen(scpname))
4780 write (iout,*) "Electrostatic potential file : ",&
4781 elename(:ilen(elename))
4782 write (iout,*) "Cumulant coefficient file : ",&
4783 fouriername(:ilen(fouriername))
4784 write (iout,*) "Torsional parameter file : ",&
4785 torname(:ilen(torname))
4786 write (iout,*) "Double torsional parameter file : ",&
4787 tordname(:ilen(tordname))
4788 write (iout,*) "SCCOR parameter file : ",&
4789 sccorname(:ilen(sccorname))
4790 write (iout,*) "Bond & inertia constant file : ",&
4791 bondname(:ilen(bondname))
4792 write (iout,*) "Bending parameter file : ",&
4793 thetname(:ilen(thetname))
4794 write (iout,*) "Rotamer parameter file : ",&
4795 rotname(:ilen(rotname))
4798 write (iout,*) "Thetpdb parameter file : ",&
4799 thetname_pdb(:ilen(thetname_pdb))
4802 write (iout,*) "Threading database : ",&
4803 patname(:ilen(patname))
4805 write (iout,*)" DIRTMP : ",&
4807 write (iout,'(80(1h-))')
4810 end subroutine openunits
4811 !-----------------------------------------------------------------------------
4814 use geometry_data, only: nres,dc
4816 ! implicit real*8 (a-h,o-z)
4817 ! include 'DIMENSIONS'
4818 ! include 'COMMON.CHAIN'
4819 ! include 'COMMON.IOUNITS'
4820 ! include 'COMMON.MD'
4823 ! real(kind=8) :: var,ene
4825 open(irest2,file=rest2name,status='unknown')
4826 read(irest2,*) totT,EK,potE,totE,t_bath
4829 ! AL 4/17/17: Now reading d_t(0,:) too
4831 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4834 ! AL 4/17/17: Now reading d_c(0,:) too
4836 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4839 read (irest2,*) iset
4843 end subroutine readrst
4844 !-----------------------------------------------------------------------------
4845 subroutine read_fragments
4849 use control_data, only:out1file
4852 ! implicit real*8 (a-h,o-z)
4853 ! include 'DIMENSIONS'
4857 ! include 'COMMON.SETUP'
4858 ! include 'COMMON.CHAIN'
4859 ! include 'COMMON.IOUNITS'
4860 ! include 'COMMON.MD'
4861 ! include 'COMMON.CONTROL'
4864 ! real(kind=8) :: var,ene
4866 read(inp,*) nset,nfrag,npair,nfrag_back
4868 !el from module energy
4869 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
4870 if(.not.allocated(wfrag_back)) then
4871 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4872 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4874 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4875 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
4877 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4878 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
4881 if(me.eq.king.or..not.out1file) &
4882 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4883 " nfrag_back",nfrag_back
4885 read(inp,*) mset(iset)
4887 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4889 if(me.eq.king.or..not.out1file) &
4890 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4891 ifrag(2,i,iset), qinfrag(i,iset)
4894 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4896 if(me.eq.king.or..not.out1file) &
4897 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4898 ipair(2,i,iset), qinpair(i,iset)
4901 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4902 wfrag_back(3,i,iset),&
4903 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4904 if(me.eq.king.or..not.out1file) &
4905 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4906 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4910 end subroutine read_fragments
4911 !-----------------------------------------------------------------------------
4913 !-----------------------------------------------------------------------------
4917 ! implicit real*8 (a-h,o-z)
4918 ! include 'DIMENSIONS'
4919 ! include 'COMMON.CSA'
4920 ! include 'COMMON.BANK'
4921 ! include 'COMMON.IOUNITS'
4923 ! integer :: ntf,ik,iw_pdb
4924 ! real(kind=8) :: var,ene
4926 open(icsa_in,file=csa_in,status="old",err=100)
4927 read(icsa_in,*) nconf
4928 read(icsa_in,*) jstart,jend
4929 read(icsa_in,*) nstmax
4930 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4931 read(icsa_in,*) nran0,nran1,irr
4932 read(icsa_in,*) nseed
4933 read(icsa_in,*) ntotal,cut1,cut2
4934 read(icsa_in,*) estop
4935 read(icsa_in,*) icmax,irestart
4936 read(icsa_in,*) ntbankm,dele,difcut
4937 read(icsa_in,*) iref,rmscut,pnccut
4938 read(icsa_in,*) ndiff
4945 end subroutine csa_read
4946 !-----------------------------------------------------------------------------
4947 subroutine initial_write
4950 ! implicit real*8 (a-h,o-z)
4951 ! include 'DIMENSIONS'
4952 ! include 'COMMON.CSA'
4953 ! include 'COMMON.BANK'
4954 ! include 'COMMON.IOUNITS'
4956 ! integer :: ntf,ik,iw_pdb
4957 ! real(kind=8) :: var,ene
4959 open(icsa_seed,file=csa_seed,status="unknown")
4960 write(icsa_seed,*) "seed"
4962 #if defined(AIX) || defined(PGI)
4963 open(icsa_history,file=csa_history,status="unknown",&
4966 open(icsa_history,file=csa_history,status="unknown",&
4969 write(icsa_history,*) nconf
4970 write(icsa_history,*) jstart,jend
4971 write(icsa_history,*) nstmax
4972 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4973 write(icsa_history,*) nran0,nran1,irr
4974 write(icsa_history,*) nseed
4975 write(icsa_history,*) ntotal,cut1,cut2
4976 write(icsa_history,*) estop
4977 write(icsa_history,*) icmax,irestart
4978 write(icsa_history,*) ntbankm,dele,difcut
4979 write(icsa_history,*) iref,rmscut,pnccut
4980 write(icsa_history,*) ndiff
4982 write(icsa_history,*)
4985 open(icsa_bank1,file=csa_bank1,status="unknown")
4986 write(icsa_bank1,*) 0
4990 end subroutine initial_write
4991 !-----------------------------------------------------------------------------
4992 subroutine restart_write
4995 ! implicit real*8 (a-h,o-z)
4996 ! include 'DIMENSIONS'
4997 ! include 'COMMON.IOUNITS'
4998 ! include 'COMMON.CSA'
4999 ! include 'COMMON.BANK'
5001 ! integer :: ntf,ik,iw_pdb
5002 ! real(kind=8) :: var,ene
5004 #if defined(AIX) || defined(PGI)
5005 open(icsa_history,file=csa_history,position="append")
5007 open(icsa_history,file=csa_history,access="append")
5009 write(icsa_history,*)
5010 write(icsa_history,*) "This is restart"
5011 write(icsa_history,*)
5012 write(icsa_history,*) nconf
5013 write(icsa_history,*) jstart,jend
5014 write(icsa_history,*) nstmax
5015 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5016 write(icsa_history,*) nran0,nran1,irr
5017 write(icsa_history,*) nseed
5018 write(icsa_history,*) ntotal,cut1,cut2
5019 write(icsa_history,*) estop
5020 write(icsa_history,*) icmax,irestart
5021 write(icsa_history,*) ntbankm,dele,difcut
5022 write(icsa_history,*) iref,rmscut,pnccut
5023 write(icsa_history,*) ndiff
5024 write(icsa_history,*)
5025 write(icsa_history,*) "irestart is: ", irestart
5027 write(icsa_history,*)
5031 end subroutine restart_write
5032 !-----------------------------------------------------------------------------
5034 !-----------------------------------------------------------------------------
5035 subroutine write_pdb(npdb,titelloc,ee)
5037 ! implicit real*8 (a-h,o-z)
5038 ! include 'DIMENSIONS'
5039 ! include 'COMMON.IOUNITS'
5040 character(len=50) :: titelloc1
5041 character*(*) :: titelloc
5042 character(len=3) :: zahl
5043 character(len=5) :: liczba5
5045 integer :: npdb !,ilen
5049 ! real(kind=8) :: var,ene
5053 if (npdb.lt.1000) then
5054 call numstr(npdb,zahl)
5055 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5057 if (npdb.lt.10000) then
5058 write(liczba5,'(i1,i4)') 0,npdb
5060 write(liczba5,'(i5)') npdb
5062 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5064 call pdbout(ee,titelloc1,ipdb)
5067 end subroutine write_pdb
5068 !-----------------------------------------------------------------------------
5070 !-----------------------------------------------------------------------------
5071 subroutine write_thread_summary
5072 ! Thread the sequence through a database of known structures
5073 use control_data, only: refstr
5075 use energy_data, only: n_ene_comp
5077 ! implicit real*8 (a-h,o-z)
5078 ! include 'DIMENSIONS'
5080 use MPI_data !include 'COMMON.INFO'
5083 ! include 'COMMON.CONTROL'
5084 ! include 'COMMON.CHAIN'
5085 ! include 'COMMON.DBASE'
5086 ! include 'COMMON.INTERACT'
5087 ! include 'COMMON.VAR'
5088 ! include 'COMMON.THREAD'
5089 ! include 'COMMON.FFIELD'
5090 ! include 'COMMON.SBRIDGE'
5091 ! include 'COMMON.HEADER'
5092 ! include 'COMMON.NAMES'
5093 ! include 'COMMON.IOUNITS'
5094 ! include 'COMMON.TIME1'
5096 integer,dimension(maxthread) :: ip
5097 real(kind=8),dimension(0:n_ene) :: energia
5099 integer :: i,j,ii,jj,ipj,ik,kk,ist
5100 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5102 write (iout,'(30x,a/)') &
5103 ' *********** Summary threading statistics ************'
5104 write (iout,'(a)') 'Initial energies:'
5105 write (iout,'(a4,2x,a12,14a14,3a8)') &
5106 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5107 'RMSnat','NatCONT','NNCONT','RMS'
5108 ! Energy sort patterns
5113 enet=ener(n_ene-1,ip(i))
5116 if (ener(n_ene-1,ip(j)).lt.enet) then
5118 enet=ener(n_ene-1,ip(j))
5130 ist=nres_base(2,ii)+ipatt(2,i)
5132 energia(i)=ener0(kk,i)
5134 etot=ener0(n_ene_comp+1,i)
5135 rmsnat=ener0(n_ene_comp+2,i)
5136 rms=ener0(n_ene_comp+3,i)
5137 frac=ener0(n_ene_comp+4,i)
5138 frac_nn=ener0(n_ene_comp+5,i)
5141 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5142 i,str_nam(ii),ist+1,&
5143 (energia(print_order(kk)),kk=1,nprint_ene),&
5144 etot,rmsnat,frac,frac_nn,rms
5146 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5147 i,str_nam(ii),ist+1,&
5148 (energia(print_order(kk)),kk=1,nprint_ene),etot
5151 write (iout,'(//a)') 'Final energies:'
5152 write (iout,'(a4,2x,a12,17a14,3a8)') &
5153 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5154 'RMSnat','NatCONT','NNCONT','RMS'
5158 ist=nres_base(2,ii)+ipatt(2,i)
5160 energia(kk)=ener(kk,ik)
5162 etot=ener(n_ene_comp+1,i)
5163 rmsnat=ener(n_ene_comp+2,i)
5164 rms=ener(n_ene_comp+3,i)
5165 frac=ener(n_ene_comp+4,i)
5166 frac_nn=ener(n_ene_comp+5,i)
5167 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5168 i,str_nam(ii),ist+1,&
5169 (energia(print_order(kk)),kk=1,nprint_ene),&
5170 etot,rmsnat,frac,frac_nn,rms
5172 write (iout,'(/a/)') 'IEXAM array:'
5173 write (iout,'(i5)') nexcl
5175 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5177 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5178 'Max. time for threading step ',max_time_for_thread,&
5179 'Average time for threading step: ',ave_time_for_thread
5181 end subroutine write_thread_summary
5182 !-----------------------------------------------------------------------------
5183 subroutine write_stat_thread(ithread,ipattern,ist)
5185 use energy_data, only: n_ene_comp
5187 ! implicit real*8 (a-h,o-z)
5188 ! include "DIMENSIONS"
5189 ! include "COMMON.CONTROL"
5190 ! include "COMMON.IOUNITS"
5191 ! include "COMMON.THREAD"
5192 ! include "COMMON.FFIELD"
5193 ! include "COMMON.DBASE"
5194 ! include "COMMON.NAMES"
5195 real(kind=8),dimension(0:n_ene) :: energia
5197 integer :: ithread,ipattern,ist,i
5198 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5200 #if defined(AIX) || defined(PGI)
5201 open(istat,file=statname,position='append')
5203 open(istat,file=statname,access='append')
5206 energia(i)=ener(i,ithread)
5208 etot=ener(n_ene_comp+1,ithread)
5209 rmsnat=ener(n_ene_comp+2,ithread)
5210 rms=ener(n_ene_comp+3,ithread)
5211 frac=ener(n_ene_comp+4,ithread)
5212 frac_nn=ener(n_ene_comp+5,ithread)
5213 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5214 ithread,str_nam(ipattern),ist+1,&
5215 (energia(print_order(i)),i=1,nprint_ene),&
5216 etot,rmsnat,frac,frac_nn,rms
5219 end subroutine write_stat_thread
5220 !-----------------------------------------------------------------------------
5222 !-----------------------------------------------------------------------------
5223 end module io_config