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)),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) :: b
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
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(ntyp)) !(ntyp)
780 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
781 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
782 allocate(msc(ntyp+1)) !(ntyp+1)
783 allocate(isc(ntyp+1)) !(ntyp+1)
784 allocate(restok(ntyp+1)) !(ntyp+1)
785 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
786 allocate(long_r_sidechain(ntyp))
787 allocate(short_r_sidechain(ntyp))
792 read (ibond,*) vbldp0,akp,mp,ip,pstok
795 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
796 dsc(i) = vbldsc0(1,i)
800 dsc_inv(i)=1.0D0/dsc(i)
804 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp,ip,pstok
806 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
807 j=1,nbondterm(i)),msc(i),isc(i),restok(i)
808 dsc(i) = vbldsc0(1,i)
812 dsc_inv(i)=1.0D0/dsc(i)
817 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
818 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
820 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
822 write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),&
823 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
825 write (iout,'(13x,3f10.5)') &
826 vbldsc0(j,i),aksc(j,i),abond0(j,i)
830 !----------------------------------------------------
831 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
832 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
833 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
834 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
835 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
836 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
846 allocate(liptranene(ntyp))
847 !C reading lipid parameters
848 write (iout,*) "iliptranpar",iliptranpar
850 read(iliptranpar,*) pepliptran
853 read(iliptranpar,*) liptranene(i)
854 print *,liptranene(i)
860 ! Read the parameters of the probability distribution/energy expression
861 ! of the virtual-bond valence angles theta
864 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
865 (bthet(j,i,1,1),j=1,2)
866 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
867 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
868 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
872 athet(1,i,1,-1)=athet(1,i,1,1)
873 athet(2,i,1,-1)=athet(2,i,1,1)
874 bthet(1,i,1,-1)=-bthet(1,i,1,1)
875 bthet(2,i,1,-1)=-bthet(2,i,1,1)
876 athet(1,i,-1,1)=-athet(1,i,1,1)
877 athet(2,i,-1,1)=-athet(2,i,1,1)
878 bthet(1,i,-1,1)=bthet(1,i,1,1)
879 bthet(2,i,-1,1)=bthet(2,i,1,1)
883 athet(1,i,-1,-1)=athet(1,-i,1,1)
884 athet(2,i,-1,-1)=-athet(2,-i,1,1)
885 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
886 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
887 athet(1,i,-1,1)=athet(1,-i,1,1)
888 athet(2,i,-1,1)=-athet(2,-i,1,1)
889 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
890 bthet(2,i,-1,1)=bthet(2,-i,1,1)
891 athet(1,i,1,-1)=-athet(1,-i,1,1)
892 athet(2,i,1,-1)=athet(2,-i,1,1)
893 bthet(1,i,1,-1)=bthet(1,-i,1,1)
894 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
899 polthet(j,i)=polthet(j,-i)
902 gthet(j,i)=gthet(j,-i)
910 'Parameters of the virtual-bond valence angles:'
911 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
912 ' ATHETA0 ',' A1 ',' A2 ',&
915 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
916 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
918 write (iout,'(/a/9x,5a/79(1h-))') &
919 'Parameters of the expression for sigma(theta_c):',&
920 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
921 ' ALPH3 ',' SIGMA0C '
923 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,&
924 (polthet(j,i),j=0,3),sigc0(i)
926 write (iout,'(/a/9x,5a/79(1h-))') &
927 'Parameters of the second gaussian:',&
928 ' THETA0 ',' SIGMA0 ',' G1 ',&
931 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),&
932 sig0(i),(gthet(j,i),j=1,3)
936 'Parameters of the virtual-bond valence angles:'
937 write (iout,'(/a/9x,5a/79(1h-))') &
938 'Coefficients of expansion',&
939 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
940 ' b1*10^1 ',' b2*10^1 '
942 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),&
943 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
944 (10*bthet(j,i,1,1),j=1,2)
946 write (iout,'(/a/9x,5a/79(1h-))') &
947 'Parameters of the expression for sigma(theta_c):',&
948 ' alpha0 ',' alph1 ',' alph2 ',&
949 ' alhp3 ',' sigma0c '
951 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),&
952 (polthet(j,i),j=0,3),sigc0(i)
954 write (iout,'(/a/9x,5a/79(1h-))') &
955 'Parameters of the second gaussian:',&
956 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
959 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),&
960 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
966 ! Read the parameters of Utheta determined from ab initio surfaces
967 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
969 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
970 ntheterm3,nsingle,ndouble
971 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
973 !----------------------------------------------------
974 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
975 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
976 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
977 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
978 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
979 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
980 !(maxtheterm,-maxthetyp1:maxthetyp1,&
981 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
982 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
983 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
984 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
985 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
986 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
987 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
988 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
989 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
990 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
991 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
992 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
993 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
994 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
995 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
996 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
997 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
999 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1001 ithetyp(i)=-ithetyp(-i)
1004 aa0thet(:,:,:,:)=0.0d0
1005 aathet(:,:,:,:,:)=0.0d0
1006 bbthet(:,:,:,:,:,:)=0.0d0
1007 ccthet(:,:,:,:,:,:)=0.0d0
1008 ddthet(:,:,:,:,:,:)=0.0d0
1009 eethet(:,:,:,:,:,:)=0.0d0
1010 ffthet(:,:,:,:,:,:,:)=0.0d0
1011 ggthet(:,:,:,:,:,:,:)=0.0d0
1013 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1015 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1016 ! VAR:1=non-glicyne non-proline 2=proline
1017 ! VAR:negative values for D-aminoacid
1019 do j=-nthetyp,nthetyp
1020 do k=-nthetyp,nthetyp
1021 read (ithep,'(6a)',end=111,err=111) res1
1022 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1023 ! VAR: aa0thet is variable describing the average value of Foureir
1024 ! VAR: expansion series
1025 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1026 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1027 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1028 read (ithep,*,end=111,err=111) &
1029 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1030 read (ithep,*,end=111,err=111) &
1031 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1032 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1033 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1034 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1036 read (ithep,*,end=111,err=111) &
1037 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1038 ffthet(lll,llll,ll,i,j,k,iblock),&
1039 ggthet(llll,lll,ll,i,j,k,iblock),&
1040 ggthet(lll,llll,ll,i,j,k,iblock),&
1041 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1046 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1047 ! coefficients of theta-and-gamma-dependent terms are zero.
1048 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1049 ! RECOMENTDED AFTER VERSION 3.3)
1053 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1054 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1056 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1057 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1060 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1062 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1065 ! AND COMMENT THE LOOPS BELOW
1069 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1070 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1072 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1073 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1076 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1078 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1083 ! Substitution for D aminoacids from symmetry.
1086 do j=-nthetyp,nthetyp
1087 do k=-nthetyp,nthetyp
1088 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1090 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1094 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1095 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1096 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1097 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1103 ffthet(llll,lll,ll,i,j,k,iblock)= &
1104 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1105 ffthet(lll,llll,ll,i,j,k,iblock)= &
1106 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1107 ggthet(llll,lll,ll,i,j,k,iblock)= &
1108 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1109 ggthet(lll,llll,ll,i,j,k,iblock)= &
1110 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1119 ! Control printout of the coefficients of virtual-bond-angle potentials
1122 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1127 write (iout,'(//4a)') &
1128 'Type ',onelett(i),onelett(j),onelett(k)
1129 write (iout,'(//a,10x,a)') " l","a[l]"
1130 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1131 write (iout,'(i2,1pe15.5)') &
1132 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1134 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1135 "b",l,"c",l,"d",l,"e",l
1137 write (iout,'(i2,4(1pe15.5))') m,&
1138 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1139 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1143 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1144 "f+",l,"f-",l,"g+",l,"g-",l
1147 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1148 ffthet(n,m,l,i,j,k,iblock),&
1149 ffthet(m,n,l,i,j,k,iblock),&
1150 ggthet(n,m,l,i,j,k,iblock),&
1151 ggthet(m,n,l,i,j,k,iblock)
1161 write (2,*) "Start reading THETA_PDB",ithep_pdb
1163 ! write (2,*) 'i=',i
1164 read (ithep_pdb,*,err=111,end=111) &
1165 a0thet(i),(athet(j,i,1,1),j=1,2),&
1166 (bthet(j,i,1,1),j=1,2)
1167 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1168 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1169 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1170 sigc0(i)=sigc0(i)**2
1173 athet(1,i,1,-1)=athet(1,i,1,1)
1174 athet(2,i,1,-1)=athet(2,i,1,1)
1175 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1176 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1177 athet(1,i,-1,1)=-athet(1,i,1,1)
1178 athet(2,i,-1,1)=-athet(2,i,1,1)
1179 bthet(1,i,-1,1)=bthet(1,i,1,1)
1180 bthet(2,i,-1,1)=bthet(2,i,1,1)
1183 a0thet(i)=a0thet(-i)
1184 athet(1,i,-1,-1)=athet(1,-i,1,1)
1185 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1186 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1187 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1188 athet(1,i,-1,1)=athet(1,-i,1,1)
1189 athet(2,i,-1,1)=-athet(2,-i,1,1)
1190 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1191 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1192 athet(1,i,1,-1)=-athet(1,-i,1,1)
1193 athet(2,i,1,-1)=athet(2,-i,1,1)
1194 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1195 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1196 theta0(i)=theta0(-i)
1200 polthet(j,i)=polthet(j,-i)
1203 gthet(j,i)=gthet(j,-i)
1206 write (2,*) "End reading THETA_PDB"
1211 !-------------------------------------------
1212 allocate(nlob(ntyp1)) !(ntyp1)
1213 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1214 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1215 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1222 gaussc(:,:,:,:)=0.0D0
1226 ! Read the parameters of the probability distribution/energy expression
1227 ! of the side chains.
1230 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1234 dsc_inv(i)=1.0D0/dsc(i)
1245 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1246 ((blower(k,l,1),l=1,k),k=1,3)
1247 censc(1,1,-i)=censc(1,1,i)
1248 censc(2,1,-i)=censc(2,1,i)
1249 censc(3,1,-i)=-censc(3,1,i)
1251 read (irotam,*,end=112,err=112) bsc(j,i)
1252 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1253 ((blower(k,l,j),l=1,k),k=1,3)
1254 censc(1,j,-i)=censc(1,j,i)
1255 censc(2,j,-i)=censc(2,j,i)
1256 censc(3,j,-i)=-censc(3,j,i)
1257 ! BSC is amplitude of Gaussian
1264 akl=akl+blower(k,m,j)*blower(l,m,j)
1268 if (((k.eq.3).and.(l.ne.3)) &
1269 .or.((l.eq.3).and.(k.ne.3))) then
1270 gaussc(k,l,j,-i)=-akl
1271 gaussc(l,k,j,-i)=-akl
1273 gaussc(k,l,j,-i)=akl
1274 gaussc(l,k,j,-i)=akl
1283 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1286 if (nlobi.gt.0) then
1288 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),&
1289 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1290 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1291 'log h',(bsc(j,i),j=1,nlobi)
1292 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1293 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1295 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1296 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1299 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1300 write (iout,'(a,f10.4,4(16x,f10.4))') &
1301 'Center ',(bsc(j,i),j=1,nlobi)
1302 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1311 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1312 ! added by Urszula Kozlowska 07/11/2007
1314 !el Maximum number of SC local term fitting function coefficiants
1315 !el integer,parameter :: maxsccoef=65
1317 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1320 read (irotam,*,end=112,err=112)
1322 read (irotam,*,end=112,err=112)
1325 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1330 ! Read the parameters of the probability distribution/energy expression
1331 ! of the side chains.
1333 write (2,*) "Start reading ROTAM_PDB"
1335 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1339 dsc_inv(i)=1.0D0/dsc(i)
1350 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1351 ((blower(k,l,1),l=1,k),k=1,3)
1353 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1354 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1355 ((blower(k,l,j),l=1,k),k=1,3)
1362 akl=akl+blower(k,m,j)*blower(l,m,j)
1372 write (2,*) "End reading ROTAM_PDB"
1378 ! Read torsional parameters in old format
1380 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1382 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1383 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1384 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1386 !el from energy module--------
1387 allocate(v1(nterm_old,ntortyp,ntortyp))
1388 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1389 !el---------------------------
1394 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1400 write (iout,'(/a/)') 'Torsional constants:'
1403 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1404 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1410 ! Read torsional parameters
1412 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1413 read (itorp,*,end=113,err=113) ntortyp
1414 !el from energy module---------
1415 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1416 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1418 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1419 allocate(vlor2(maxlor,ntortyp,ntortyp))
1420 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1421 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1423 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1424 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1425 !el---------------------------
1428 !el---------------------------
1430 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1432 itortyp(i)=-itortyp(-i)
1434 itortyp(ntyp1)=ntortyp
1435 itortyp(-ntyp1)=-ntortyp
1437 write (iout,*) 'ntortyp',ntortyp
1439 do j=-ntortyp+1,ntortyp-1
1440 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1442 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1443 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1446 do k=1,nterm(i,j,iblock)
1447 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1449 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1450 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1451 v0ij=v0ij+si*v1(k,i,j,iblock)
1453 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1454 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1455 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1457 do k=1,nlor(i,j,iblock)
1458 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1459 vlor2(k,i,j),vlor3(k,i,j)
1460 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1463 v0(-i,-j,iblock)=v0ij
1469 write (iout,'(/a/)') 'Torsional constants:'
1471 do i=-ntortyp,ntortyp
1472 do j=-ntortyp,ntortyp
1473 write (iout,*) 'ityp',i,' jtyp',j
1474 write (iout,*) 'Fourier constants'
1475 do k=1,nterm(i,j,iblock)
1476 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1479 write (iout,*) 'Lorenz constants'
1480 do k=1,nlor(i,j,iblock)
1481 write (iout,'(3(1pe15.5))') &
1482 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1488 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1490 ! 6/23/01 Read parameters for double torsionals
1492 !el from energy module------------
1493 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1494 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1495 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1496 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1497 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1498 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1499 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1500 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1501 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1502 !---------------------------------
1506 do j=-ntortyp+1,ntortyp-1
1507 do k=-ntortyp+1,ntortyp-1
1508 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1509 ! write (iout,*) "OK onelett",
1512 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1513 .or. t3.ne.toronelet(k)) then
1514 write (iout,*) "Error in double torsional parameter file",&
1517 call MPI_Finalize(Ierror)
1519 stop "Error in double torsional parameter file"
1521 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1522 ntermd_2(i,j,k,iblock)
1523 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1524 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1525 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1526 ntermd_1(i,j,k,iblock))
1527 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1528 ntermd_1(i,j,k,iblock))
1529 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1530 ntermd_1(i,j,k,iblock))
1531 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1532 ntermd_1(i,j,k,iblock))
1533 ! Martix of D parameters for one dimesional foureir series
1534 do l=1,ntermd_1(i,j,k,iblock)
1535 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1536 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1537 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1538 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1539 ! write(iout,*) "whcodze" ,
1540 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1542 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1543 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1544 v2s(m,l,i,j,k,iblock),&
1545 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1546 ! Martix of D parameters for two dimesional fourier series
1547 do l=1,ntermd_2(i,j,k,iblock)
1549 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1550 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1551 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1552 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1561 write (iout,*) 'Constants for double torsionals'
1564 do j=-ntortyp+1,ntortyp-1
1565 do k=-ntortyp+1,ntortyp-1
1566 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1567 ' nsingle',ntermd_1(i,j,k,iblock),&
1568 ' ndouble',ntermd_2(i,j,k,iblock)
1570 write (iout,*) 'Single angles:'
1571 do l=1,ntermd_1(i,j,k,iblock)
1572 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1573 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1574 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1575 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1578 write (iout,*) 'Pairs of angles:'
1579 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1580 do l=1,ntermd_2(i,j,k,iblock)
1581 write (iout,'(i5,20f10.5)') &
1582 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1585 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1586 do l=1,ntermd_2(i,j,k,iblock)
1587 write (iout,'(i5,20f10.5)') &
1588 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1589 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1598 ! Read of Side-chain backbone correlation parameters
1599 ! Modified 11 May 2012 by Adasko
1602 read (isccor,*,end=119,err=119) nsccortyp
1604 !el from module energy-------------
1605 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1606 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1607 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1608 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1609 !-----------------------------------
1611 !el from module energy-------------
1612 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1614 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1616 isccortyp(i)=-isccortyp(-i)
1618 iscprol=isccortyp(20)
1619 ! write (iout,*) 'ntortyp',ntortyp
1621 !c maxinter is maximum interaction sites
1622 !el from module energy---------
1623 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1624 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1625 -nsccortyp:nsccortyp))
1626 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1627 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1628 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1629 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1630 !-----------------------------------
1634 read (isccor,*,end=119,err=119) &
1635 nterm_sccor(i,j),nlor_sccor(i,j)
1641 nterm_sccor(-i,j)=nterm_sccor(i,j)
1642 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1643 nterm_sccor(i,-j)=nterm_sccor(i,j)
1644 do k=1,nterm_sccor(i,j)
1645 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1647 if (j.eq.iscprol) then
1648 if (i.eq.isccortyp(10)) then
1649 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1650 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1652 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1653 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1654 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1655 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1656 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1657 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1658 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1659 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1662 if (i.eq.isccortyp(10)) then
1663 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1664 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1666 if (j.eq.isccortyp(10)) then
1667 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1668 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1670 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1671 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1672 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1673 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1674 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1675 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1679 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1680 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1681 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1682 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1685 do k=1,nlor_sccor(i,j)
1686 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1687 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1688 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1689 (1+vlor3sccor(k,i,j)**2)
1691 v0sccor(l,i,j)=v0ijsccor
1692 v0sccor(l,-i,j)=v0ijsccor1
1693 v0sccor(l,i,-j)=v0ijsccor2
1694 v0sccor(l,-i,-j)=v0ijsccor3
1700 !el from module energy-------------
1701 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1703 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1704 ! write (iout,*) 'ntortyp',ntortyp
1706 !c maxinter is maximum interaction sites
1707 !el from module energy---------
1708 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1709 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1710 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1711 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1712 !-----------------------------------
1716 read (isccor,*,end=119,err=119) &
1717 nterm_sccor(i,j),nlor_sccor(i,j)
1721 do k=1,nterm_sccor(i,j)
1722 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1724 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1727 do k=1,nlor_sccor(i,j)
1728 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1729 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1730 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1731 (1+vlor3sccor(k,i,j)**2)
1733 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1741 write (iout,'(/a/)') 'Torsional constants:'
1744 write (iout,*) 'ityp',i,' jtyp',j
1745 write (iout,*) 'Fourier constants'
1746 do k=1,nterm_sccor(i,j)
1747 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1749 write (iout,*) 'Lorenz constants'
1750 do k=1,nlor_sccor(i,j)
1751 write (iout,'(3(1pe15.5))') &
1752 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1759 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1760 ! interaction energy of the Gly, Ala, and Pro prototypes.
1765 write (iout,*) "Coefficients of the cumulants"
1767 read (ifourier,*) nloctyp
1768 !write(iout,*) "nloctyp",nloctyp
1769 !el from module energy-------
1770 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1771 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1772 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1773 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1774 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1775 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1776 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1777 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1781 !--------------------------------
1784 read (ifourier,*,end=115,err=115)
1785 read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
1787 write (iout,*) 'Type',i
1788 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
1798 B1tilde(1,-i) =-b(3)
1800 ! b1tilde(1,i)=0.0d0
1801 ! b1tilde(2,i)=0.0d0
1826 Ctilde(1,2,-i)=-b(9)
1830 ! Ctilde(1,1,i)=0.0d0
1831 ! Ctilde(1,2,i)=0.0d0
1832 ! Ctilde(2,1,i)=0.0d0
1833 ! Ctilde(2,2,i)=0.0d0
1851 Dtilde(1,2,-i)=-b(8)
1855 ! Dtilde(1,1,i)=0.0d0
1856 ! Dtilde(1,2,i)=0.0d0
1857 ! Dtilde(2,1,i)=0.0d0
1858 ! Dtilde(2,2,i)=0.0d0
1859 EE(1,1,i)= b(10)+b(11)
1860 EE(2,2,i)=-b(10)+b(11)
1861 EE(2,1,i)= b(12)-b(13)
1862 EE(1,2,i)= b(12)+b(13)
1863 EE(1,1,-i)= b(10)+b(11)
1864 EE(2,2,-i)=-b(10)+b(11)
1865 EE(2,1,-i)=-b(12)+b(13)
1866 EE(1,2,-i)=-b(12)-b(13)
1872 ! ee(2,1,i)=ee(1,2,i)
1876 write (iout,*) 'Type',i
1878 write(iout,*) B1(1,i),B1(2,i)
1880 write(iout,*) B2(1,i),B2(2,i)
1883 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
1887 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
1891 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
1896 ! Read electrostatic-interaction parameters
1901 write (iout,'(/a)') 'Electrostatic interaction constants:'
1902 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
1903 'IT','JT','APP','BPP','AEL6','AEL3'
1905 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
1906 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
1907 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
1908 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
1913 app (i,j)=epp(i,j)*rri*rri
1914 bpp (i,j)=-2.0D0*epp(i,j)*rri
1915 ael6(i,j)=elpp6(i,j)*4.2D0**6
1916 ael3(i,j)=elpp3(i,j)*4.2D0**3
1918 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
1924 ! Read side-chain interaction parameters.
1926 !el from module energy - COMMON.INTERACT-------
1927 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
1928 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
1929 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
1930 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
1931 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
1932 allocate(epslip(ntyp,ntyp))
1940 !--------------------------------
1942 read (isidep,*,end=117,err=117) ipot,expon
1943 if (ipot.lt.1 .or. ipot.gt.5) then
1944 write (iout,'(2a)') 'Error while reading SC interaction',&
1945 'potential file - unknown potential type.'
1947 call MPI_Finalize(Ierror)
1953 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
1954 ', exponents are ',expon,2*expon
1955 ! goto (10,20,30,30,40) ipot
1957 !----------------------- LJ potential ---------------------------------
1959 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1960 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1961 (sigma0(i),i=1,ntyp)
1963 write (iout,'(/a/)') 'Parameters of the LJ potential:'
1964 write (iout,'(a/)') 'The epsilon array:'
1965 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1966 write (iout,'(/a)') 'One-body parameters:'
1967 write (iout,'(a,4x,a)') 'residue','sigma'
1968 write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
1971 !----------------------- LJK potential --------------------------------
1973 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1974 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
1975 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
1977 write (iout,'(/a/)') 'Parameters of the LJK potential:'
1978 write (iout,'(a/)') 'The epsilon array:'
1979 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
1980 write (iout,'(/a)') 'One-body parameters:'
1981 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
1982 write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
1986 !---------------------- GB or BP potential -----------------------------
1990 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
1992 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
1993 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
1994 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
1995 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
1997 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2000 ! For the GB potential convert sigma'**2 into chi'
2003 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2007 write (iout,'(/a/)') 'Parameters of the BP potential:'
2008 write (iout,'(a/)') 'The epsilon array:'
2009 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2010 write (iout,'(/a)') 'One-body parameters:'
2011 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2013 write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
2014 chip(i),alp(i),i=1,ntyp)
2017 !--------------------- GBV potential -----------------------------------
2019 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2020 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2021 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2022 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2024 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2025 write (iout,'(a/)') 'The epsilon array:'
2026 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2027 write (iout,'(/a)') 'One-body parameters:'
2028 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2029 's||/s_|_^2',' chip ',' alph '
2030 write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
2031 sigii(i),chip(i),alp(i),i=1,ntyp)
2034 write(iout,*)"Wrong ipot"
2039 !-----------------------------------------------------------------------
2040 ! Calculate the "working" parameters of SC interactions.
2042 !el from module energy - COMMON.INTERACT-------
2043 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2044 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2045 allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2046 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2048 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2061 sc_aa_tube_par(:)=0.0d0
2062 sc_bb_tube_par(:)=0.0d0
2064 !--------------------------------
2069 epslip(i,j)=epslip(j,i)
2074 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2075 sigma(j,i)=sigma(i,j)
2076 rs0(i,j)=dwa16*sigma(i,j)
2080 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2081 'Working parameters of the SC interactions:',&
2082 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2087 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2096 sigeps=dsign(1.0D0,epsij)
2098 aa_aq(i,j)=epsij*rrij*rrij
2099 bb_aq(i,j)=-sigeps*epsij*rrij
2100 aa_aq(j,i)=aa_aq(i,j)
2101 bb_aq(j,i)=bb_aq(i,j)
2102 epsijlip=epslip(i,j)
2103 sigeps=dsign(1.0D0,epsijlip)
2104 epsijlip=dabs(epsijlip)
2105 aa_lip(i,j)=epsijlip*rrij*rrij
2106 bb_lip(i,j)=-sigeps*epsijlip*rrij
2107 aa_lip(j,i)=aa_lip(i,j)
2108 bb_lip(j,i)=bb_lip(i,j)
2109 !C write(iout,*) aa_lip
2111 sigt1sq=sigma0(i)**2
2112 sigt2sq=sigma0(j)**2
2115 ratsig1=sigt2sq/sigt1sq
2116 ratsig2=1.0D0/ratsig1
2117 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2118 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2119 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2123 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2124 sigmaii(i,j)=rsum_max
2125 sigmaii(j,i)=rsum_max
2127 ! sigmaii(i,j)=r0(i,j)
2128 ! sigmaii(j,i)=r0(i,j)
2130 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2131 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2132 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2133 augm(i,j)=epsij*r_augm**(2*expon)
2134 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2141 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2142 restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2143 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2147 write(iout,*) "tube param"
2148 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2149 ccavtubpep,dcavtubpep,tubetranenepep
2150 sigmapeptube=sigmapeptube**6
2151 sigeps=dsign(1.0D0,epspeptube)
2152 epspeptube=dabs(epspeptube)
2153 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2154 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2155 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2157 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2158 ccavtub(i),dcavtub(i),tubetranene(i)
2159 sigmasctube=sigmasctube**6
2160 sigeps=dsign(1.0D0,epssctube)
2161 epssctube=dabs(epssctube)
2162 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2163 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2164 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2167 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2172 ! Define the SC-p interaction constants (hard-coded; old style)
2175 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2177 ! aad(i,1)=0.3D0*4.0D0**12
2178 ! Following line for constants currently implemented
2179 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2180 aad(i,1)=1.5D0*4.0D0**12
2181 ! aad(i,1)=0.17D0*5.6D0**12
2183 ! "Soft" SC-p repulsion
2185 ! Following line for constants currently implemented
2186 ! aad(i,1)=0.3D0*4.0D0**6
2187 ! "Hard" SC-p repulsion
2188 bad(i,1)=3.0D0*4.0D0**6
2189 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2198 ! 8/9/01 Read the SC-p interaction constants from file
2201 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2204 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2205 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2206 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2207 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2211 write (iout,*) "Parameters of SC-p interactions:"
2213 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2214 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2220 ! Define the constants of the disulfide bridge
2224 ! Old arbitrary potential - commented out.
2229 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2230 ! energy surface of diethyl disulfide.
2231 ! A. Liwo and U. Kozlowska, 11/24/03
2248 write (iout,'(/a)') "Disulfide bridge parameters:"
2249 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2250 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2251 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2252 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2256 111 write (iout,*) "Error reading bending energy parameters."
2258 112 write (iout,*) "Error reading rotamer energy parameters."
2260 113 write (iout,*) "Error reading torsional energy parameters."
2262 114 write (iout,*) "Error reading double torsional energy parameters."
2264 115 write (iout,*) &
2265 "Error reading cumulant (multibody energy) parameters."
2267 116 write (iout,*) "Error reading electrostatic energy parameters."
2269 117 write (iout,*) "Error reading side chain interaction parameters."
2271 118 write (iout,*) "Error reading SCp interaction parameters."
2273 119 write (iout,*) "Error reading SCCOR parameters"
2276 call MPI_Finalize(Ierror)
2280 end subroutine parmread
2282 !-----------------------------------------------------------------------------
2284 !-----------------------------------------------------------------------------
2285 subroutine printmat(ldim,m,n,iout,key,a)
2288 character(len=3),dimension(n) :: key
2289 real(kind=8),dimension(ldim,n) :: a
2291 integer :: i,j,k,m,iout,nlim
2295 write (iout,1000) (key(k),k=i,nlim)
2297 1000 format (/5x,8(6x,a3))
2298 1020 format (/80(1h-)/)
2300 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2303 1010 format (a3,2x,8(f9.4))
2305 end subroutine printmat
2306 !-----------------------------------------------------------------------------
2308 !-----------------------------------------------------------------------------
2310 ! Read the PDB file and convert the peptide geometry into virtual-chain
2313 use energy_data, only: itype
2317 use control, only: rescode
2318 ! implicit real*8 (a-h,o-z)
2319 ! include 'DIMENSIONS'
2320 ! include 'COMMON.LOCAL'
2321 ! include 'COMMON.VAR'
2322 ! include 'COMMON.CHAIN'
2323 ! include 'COMMON.INTERACT'
2324 ! include 'COMMON.IOUNITS'
2325 ! include 'COMMON.GEO'
2326 ! include 'COMMON.NAMES'
2327 ! include 'COMMON.CONTROL'
2328 ! include 'COMMON.DISTFIT'
2329 ! include 'COMMON.SETUP'
2330 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,&
2332 logical :: lprn=.true.,fail
2333 real(kind=8),dimension(3) :: e1,e2,e3
2334 real(kind=8) :: dcj,efree_temp
2335 character(len=3) :: seq,res
2336 character(len=5) :: atom
2337 character(len=80) :: card
2338 real(kind=8),dimension(3,20) :: sccor
2339 integer :: kkk,lll,icha,kupa !rescode,
2342 integer,dimension(2,maxres/3) :: hfrag_alloc
2343 integer,dimension(4,maxres/3) :: bfrag_alloc
2344 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2350 ! write (2,*) "UNRES_PDB",unres_pdb
2358 !-----------------------------
2359 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2360 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2363 read (ipdbin,'(a80)',end=10) card
2364 ! write (iout,'(a)') card
2365 if (card(:5).eq.'HELIX') then
2368 read(card(22:25),*) hfrag(1,nhfrag)
2369 read(card(34:37),*) hfrag(2,nhfrag)
2371 if (card(:5).eq.'SHEET') then
2374 read(card(24:26),*) bfrag(1,nbfrag)
2375 read(card(35:37),*) bfrag(2,nbfrag)
2376 !rc----------------------------------------
2377 !rc to be corrected !!!
2378 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2379 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2380 !rc----------------------------------------
2382 if (card(:3).eq.'END') then
2384 else if (card(:3).eq.'TER') then
2388 itype(ires_old)=ntyp1
2389 itype(ires_old-1)=ntyp1
2391 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2394 dc(j,ires)=sccor(j,iii)
2397 call sccenter(ires,iii,sccor)
2403 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2404 ! Fish out the ATOM cards.
2405 if (index(card(1:4),'ATOM').gt.0) then
2406 read (card(12:16),*) atom
2407 ! write (iout,*) "! ",atom," !",ires
2408 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2409 read (card(23:26),*) ires
2410 read (card(18:20),'(a3)') res
2411 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2412 ! & " ires_old",ires_old
2413 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2414 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2415 if (ires-ishift+ishift1.ne.ires_old) then
2416 ! Calculate the CM of the preceding residue.
2417 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2419 ! write (iout,*) "Calculating sidechain center iii",iii
2422 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2425 call sccenter(ires_old,iii,sccor)
2429 ! Start new residue.
2430 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2433 else if (ibeg.eq.1) then
2434 write (iout,*) "BEG ires",ires
2436 if (res.ne.'GLY' .and. res.ne. 'ACE') then
2440 ires=ires-ishift+ishift1
2442 ! write (iout,*) "ishift",ishift," ires",ires,&
2443 ! " ires_old",ires_old
2445 else if (ibeg.eq.2) then
2447 ishift=-ires_old+ires-1 !!!!!
2448 ishift1=ishift1-1 !!!!!
2449 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2450 ires=ires-ishift+ishift1
2454 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2455 ires=ires-ishift+ishift1
2458 if (res.eq.'ACE' .or. res.eq.'NHE') then
2461 itype(ires)=rescode(ires,res,0)
2464 ires=ires-ishift+ishift1
2466 ! write (iout,*) "ires_old",ires_old," ires",ires
2467 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
2470 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
2471 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
2472 res.eq.'NHE'.and.atom(:2).eq.'HN') then
2473 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
2474 ! write (iout,*) "backbone ",atom
2476 write (iout,'(2i3,2x,a,3f8.3)') &
2477 ires,itype(ires),res,(c(j,ires),j=1,3)
2481 sccor(j,iii)=c(j,ires)
2483 ! write (*,*) card(23:27),ires,itype(ires)
2484 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
2485 atom.ne.'N' .and. atom.ne.'C' .and. &
2486 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
2487 atom.ne.'OXT' .and. atom(:2).ne.'3H') then
2488 ! write (iout,*) "sidechain ",atom
2490 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
2494 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
2495 if (ires.eq.0) return
2496 ! Calculate dummy residue coordinates inside the "chain" of a multichain
2500 ! write (iout,*) i,itype(i)
2501 ! if (itype(i).eq.ntyp1) then
2502 ! write (iout,*) "dummy",i,itype(i)
2504 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
2505 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
2509 if (itype(i).eq.ntyp1) then
2510 if (itype(i+1).eq.ntyp1) then
2511 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
2512 ! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
2513 ! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
2515 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2516 ! print *,i,'tu dochodze'
2517 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
2525 c(j,i)=c(j,i-1)-1.9d0*e2(j)
2529 dcj=(c(j,i-2)-c(j,i-3))/2.0
2530 if (dcj.eq.0) dcj=1.23591524223
2535 else !itype(i+1).eq.ntyp1
2537 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2538 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
2545 c(j,i)=c(j,i+1)-1.9d0*e2(j)
2549 dcj=(c(j,i+3)-c(j,i+2))/2.0
2550 if (dcj.eq.0) dcj=1.23591524223
2555 endif !itype(i+1).eq.ntyp1
2556 endif !itype.eq.ntyp1
2559 ! Calculate the CM of the last side chain.
2563 dc(j,ires)=sccor(j,iii)
2566 call sccenter(ires,iii,sccor)
2572 if (itype(nres).ne.10) then
2576 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
2577 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
2584 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
2588 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
2589 c(j,nres)=c(j,nres-1)+dcj
2590 c(j,2*nres)=c(j,nres)
2594 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
2596 if (nres.ne.nres0) then
2597 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
2599 stop "Error nres value in WHAM input"
2602 !---------------------------------
2603 !el reallocate tables
2606 ! hfrag_alloc(j,i)=hfrag(j,i)
2609 ! bfrag_alloc(j,i)=bfrag(j,i)
2615 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
2616 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
2617 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
2618 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
2622 ! hfrag(j,i)=hfrag_alloc(j,i)
2627 ! bfrag(j,i)=bfrag_alloc(j,i)
2630 !el end reallocate tables
2631 !---------------------------------
2639 c(j,2*nres)=c(j,nres)
2641 if (itype(1).eq.ntyp1) then
2645 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
2646 call refsys(2,3,4,e1,e2,e3,fail)
2653 c(j,1)=c(j,2)-1.9d0*e2(j)
2657 dcj=(c(j,4)-c(j,3))/2.0
2663 ! Copy the coordinates to reference coordinates
2669 ! Calculate internal coordinates.
2671 write (iout,'(/a)') &
2672 "Cartesian coordinates of the reference structure"
2673 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
2674 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
2676 write (iout,'(a3,1x,i3,3f8.3,5x,3f8.3)') &
2677 restyp(itype(ires)),ires,(c(j,ires),j=1,3),&
2678 (c(j,ires+nres),j=1,3)
2681 ! znamy już nres wiec mozna alokowac tablice
2682 ! Calculate internal coordinates.
2683 if(me.eq.king.or..not.out1file)then
2684 write (iout,'(a)') &
2685 "Backbone and SC coordinates as read from the PDB"
2687 write (iout,'(2i3,2x,a,3f8.3,5x,3f8.3)') &
2688 ires,itype(ires),restyp(itype(ires)),(c(j,ires),j=1,3),&
2689 (c(j,nres+ires),j=1,3)
2693 if(.not.allocated(vbld)) then
2694 allocate(vbld(2*nres))
2699 if(.not.allocated(vbld_inv)) then
2700 allocate(vbld_inv(2*nres))
2706 if(.not.allocated(theta)) then
2707 allocate(theta(nres+2))
2711 if(.not.allocated(phi)) allocate(phi(nres+2))
2712 if(.not.allocated(alph)) allocate(alph(nres+2))
2713 if(.not.allocated(omeg)) allocate(omeg(nres+2))
2714 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
2715 if(.not.allocated(phiref)) allocate(phiref(nres+2))
2716 if(.not.allocated(costtab)) allocate(costtab(nres))
2717 if(.not.allocated(sinttab)) allocate(sinttab(nres))
2718 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
2719 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
2720 if(.not.allocated(xxref)) allocate(xxref(nres))
2721 if(.not.allocated(yyref)) allocate(yyref(nres))
2722 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
2723 if(.not.allocated(dc_norm)) then
2724 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
2725 allocate(dc_norm(3,0:2*nres+2))
2729 call int_from_cart(.true.,.false.)
2730 call sc_loc_geom(.false.)
2732 thetaref(i)=theta(i)
2742 dc(j,i)=c(j,i+1)-c(j,i)
2743 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
2748 dc(j,i+nres)=c(j,i+nres)-c(j,i)
2749 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
2751 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
2755 ! Copy the coordinates to reference coordinates
2756 ! Splits to single chain if occurs
2760 ! cref(j,i,cou)=c(j,i)
2764 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
2765 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
2766 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
2767 !-----------------------------
2773 !c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2775 if ((itype(i-1).eq.ntyp1).and.(i.gt.2)) then
2778 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
2783 cref(j,i,cou)=c(j,i)
2784 cref(j,i+nres,cou)=c(j,i+nres)
2786 chain_rep(j,lll,kkk)=c(j,i)
2787 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
2791 write (iout,*) chain_length
2792 if (chain_length.eq.0) chain_length=nres
2794 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
2795 chain_rep(j,chain_length+nres,symetr) &
2796 =chain_rep(j,chain_length+nres,1)
2799 ! write (iout,*) "spraw lancuchy",chain_length,symetr
2801 ! do kkk=1,chain_length
2802 ! write (iout,*) itype(kkk),(chain_rep(j,kkk,i), j=1,3)
2806 ! makes copy of chains
2807 write (iout,*) "symetr", symetr
2812 if (symetr.gt.1) then
2819 write(iout,*) (tabperm(i,kkk),kkk=1,4)
2825 ! write (iout,*) i,icha
2826 do lll=1,chain_length
2828 if (cou.le.nres) then
2830 kupa=mod(lll,chain_length)
2831 iprzes=(kkk-1)*chain_length+lll
2832 if (kupa.eq.0) kupa=chain_length
2833 ! write (iout,*) "kupa", kupa
2834 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
2835 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
2842 !-koniec robienia kopii
2845 write (iout,*) "nowa struktura", nperm
2847 write (iout,110) restyp(itype(i)),i,cref(1,i,kkk),&
2849 cref(3,i,kkk),cref(1,nres+i,kkk),&
2850 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
2852 100 format (//' alpha-carbon coordinates ',&
2853 ' centroid coordinates'/ &
2854 ' ', 6X,'X',11X,'Y',11X,'Z', &
2855 10X,'X',11X,'Y',11X,'Z')
2856 110 format (a,'(',i3,')',6f12.5)
2862 bfrag(i,j)=bfrag(i,j)-ishift
2868 hfrag(i,j)=hfrag(i,j)-ishift
2874 end subroutine readpdb
2875 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
2876 !-----------------------------------------------------------------------------
2878 !-----------------------------------------------------------------------------
2879 subroutine read_control
2893 use random, only: random_init
2894 ! implicit real*8 (a-h,o-z)
2895 ! include 'DIMENSIONS'
2897 use prng, only:prng_restart
2899 logical :: OKRandom!, prng_restart
2902 ! include 'COMMON.IOUNITS'
2903 ! include 'COMMON.TIME1'
2904 ! include 'COMMON.THREAD'
2905 ! include 'COMMON.SBRIDGE'
2906 ! include 'COMMON.CONTROL'
2907 ! include 'COMMON.MCM'
2908 ! include 'COMMON.MAP'
2909 ! include 'COMMON.HEADER'
2910 ! include 'COMMON.CSA'
2911 ! include 'COMMON.CHAIN'
2912 ! include 'COMMON.MUCA'
2913 ! include 'COMMON.MD'
2914 ! include 'COMMON.FFIELD'
2915 ! include 'COMMON.INTERACT'
2916 ! include 'COMMON.SETUP'
2917 !el integer :: KDIAG,ICORFL,IXDR
2918 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
2919 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
2920 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
2921 ! character(len=80) :: ucase
2922 character(len=640) :: controlcard
2924 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
2930 read (INP,'(a)') titel
2931 call card_concat(controlcard,.true.)
2932 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
2933 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
2934 call reada(controlcard,'SEED',seed,0.0D0)
2935 call random_init(seed)
2936 ! Set up the time limit (caution! The time must be input in minutes!)
2937 read_cart=index(controlcard,'READ_CART').gt.0
2938 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
2939 call readi(controlcard,'SYM',symetr,1)
2940 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
2941 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
2942 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
2943 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
2944 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
2945 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
2946 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
2947 call reada(controlcard,'DRMS',drms,0.1D0)
2948 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2949 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
2950 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
2951 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
2952 write (iout,'(a,f10.1)')'DRMS = ',drms
2953 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
2954 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
2956 call readi(controlcard,'NZ_START',nz_start,0)
2957 call readi(controlcard,'NZ_END',nz_end,0)
2958 call readi(controlcard,'IZ_SC',iz_sc,0)
2959 timlim=60.0D0*timlim
2960 safety = 60.0d0*safety
2963 call reada(controlcard,"T_BATH",t_bath,300.0d0)
2964 !C SHIELD keyword sets if the shielding effect of side-chains is used
2965 !C 0 denots no shielding is used all peptide are equally despite the
2966 !C solvent accesible area
2967 !C 1 the newly introduced function
2968 !C 2 reseved for further possible developement
2969 call readi(controlcard,'SHIELD',shield_mode,0)
2970 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2971 write(iout,*) "shield_mode",shield_mode
2972 !C Varibles set size of box
2973 call reada(controlcard,'BOXX',boxxsize,100.0d0)
2974 call reada(controlcard,'BOXY',boxysize,100.0d0)
2975 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
2976 call readi(controlcard,'TUBEMOD',tubemode,0)
2977 write (iout,*) TUBEmode,"TUBEMODE"
2978 if (TUBEmode.gt.0) then
2979 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
2980 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
2981 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
2982 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
2983 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
2984 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
2985 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
2986 buftubebot=bordtubebot+tubebufthick
2987 buftubetop=bordtubetop-tubebufthick
2990 ! CUTOFFF ON ELECTROSTATICS
2991 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
2992 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
2993 write(iout,*) "R_CUT_ELE=",r_cut_ele
2994 ! Lipidic parameters
2995 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
2996 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
2997 if (lipthick.gt.0.0d0) then
2998 bordliptop=(boxzsize+lipthick)/2.0
2999 bordlipbot=bordliptop-lipthick
3000 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3001 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3002 buflipbot=bordlipbot+lipbufthick
3003 bufliptop=bordliptop-lipbufthick
3004 if ((lipbufthick*2.0d0).gt.lipthick) &
3005 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3006 endif !lipthick.gt.0
3007 write(iout,*) "bordliptop=",bordliptop
3008 write(iout,*) "bordlipbot=",bordlipbot
3009 write(iout,*) "bufliptop=",bufliptop
3010 write(iout,*) "buflipbot=",buflipbot
3011 write (iout,*) "SHIELD MODE",shield_mode
3013 !C-------------------------
3014 minim=(index(controlcard,'MINIMIZE').gt.0)
3015 dccart=(index(controlcard,'CART').gt.0)
3016 overlapsc=(index(controlcard,'OVERLAP').gt.0)
3017 overlapsc=.not.overlapsc
3018 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3019 searchsc=.not.searchsc
3020 sideadd=(index(controlcard,'SIDEADD').gt.0)
3021 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3022 outpdb=(index(controlcard,'PDBOUT').gt.0)
3023 outmol2=(index(controlcard,'MOL2OUT').gt.0)
3024 pdbref=(index(controlcard,'PDBREF').gt.0)
3025 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3026 indpdb=index(controlcard,'PDBSTART')
3027 extconf=(index(controlcard,'EXTCONF').gt.0)
3028 call readi(controlcard,'IPRINT',iprint,0)
3029 call readi(controlcard,'MAXGEN',maxgen,10000)
3030 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3031 call readi(controlcard,"KDIAG",kdiag,0)
3032 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3033 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3034 write (iout,*) "RESCALE_MODE",rescale_mode
3035 split_ene=index(controlcard,'SPLIT_ENE').gt.0
3036 if (index(controlcard,'REGULAR').gt.0.0D0) then
3037 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3041 if (index(controlcard,'CHECKGRAD').gt.0) then
3043 if (index(controlcard,'CART').gt.0) then
3045 elseif (index(controlcard,'CARINT').gt.0) then
3050 elseif (index(controlcard,'THREAD').gt.0) then
3052 call readi(controlcard,'THREAD',nthread,0)
3053 if (nthread.gt.0) then
3054 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3057 write (iout,'(a)')'A number has to follow the THREAD keyword.'
3058 stop 'Error termination in Read_Control.'
3060 else if (index(controlcard,'MCMA').gt.0) then
3062 else if (index(controlcard,'MCEE').gt.0) then
3064 else if (index(controlcard,'MULTCONF').gt.0) then
3066 else if (index(controlcard,'MAP').gt.0) then
3068 call readi(controlcard,'MAP',nmap,0)
3069 else if (index(controlcard,'CSA').gt.0) then
3071 !rc else if (index(controlcard,'ZSCORE').gt.0) then
3073 !rc ZSCORE is rm from UNRES, modecalc=9 is available
3076 !fcm else if (index(controlcard,'MCMF').gt.0) then
3078 else if (index(controlcard,'SOFTREG').gt.0) then
3080 else if (index(controlcard,'CHECK_BOND').gt.0) then
3082 else if (index(controlcard,'TEST').gt.0) then
3084 else if (index(controlcard,'MD').gt.0) then
3086 else if (index(controlcard,'RE ').gt.0) then
3090 lmuca=index(controlcard,'MUCA').gt.0
3091 call readi(controlcard,'MUCADYN',mucadyn,0)
3092 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3093 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3095 write (iout,*) 'MUCADYN=',mucadyn
3096 write (iout,*) 'MUCASMOOTH=',muca_smooth
3099 iscode=index(controlcard,'ONE_LETTER')
3100 indphi=index(controlcard,'PHI')
3101 indback=index(controlcard,'BACK')
3102 iranconf=index(controlcard,'RAND_CONF')
3103 i2ndstr=index(controlcard,'USE_SEC_PRED')
3104 gradout=index(controlcard,'GRADOUT').gt.0
3105 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3106 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3107 if (me.eq.king .or. .not.out1file ) &
3108 write (iout,*) "DISTCHAINMAX",distchainmax
3110 if(me.eq.king.or..not.out1file) &
3111 write (iout,'(2a)') diagmeth(kdiag),&
3112 ' routine used to diagonalize matrices.'
3113 if (shield_mode.gt.0) then
3115 !C VSolvSphere the volume of solving sphere
3117 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3118 !C there will be no distinction between proline peptide group and normal peptide
3119 !C group in case of shielding parameters
3120 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3121 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3122 write (iout,*) VSolvSphere,VSolvSphere_div
3123 !C long axis of side chain
3125 long_r_sidechain(i)=vbldsc0(1,i)
3126 short_r_sidechain(i)=sigma0(i)
3127 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3133 end subroutine read_control
3134 !-----------------------------------------------------------------------------
3135 subroutine read_REMDpar
3137 ! Read REMD settings
3144 use control_data, only:out1file
3145 ! implicit real*8 (a-h,o-z)
3146 ! include 'DIMENSIONS'
3147 ! include 'COMMON.IOUNITS'
3148 ! include 'COMMON.TIME1'
3149 ! include 'COMMON.MD'
3152 !el include 'COMMON.LANGEVIN'
3154 !el include 'COMMON.LANGEVIN.lang0'
3156 ! include 'COMMON.INTERACT'
3157 ! include 'COMMON.NAMES'
3158 ! include 'COMMON.GEO'
3159 ! include 'COMMON.REMD'
3160 ! include 'COMMON.CONTROL'
3161 ! include 'COMMON.SETUP'
3162 ! character(len=80) :: ucase
3163 character(len=320) :: controlcard
3164 character(len=3200) :: controlcard1
3165 integer :: iremd_m_total
3168 ! real(kind=8) :: var,ene
3170 if(me.eq.king.or..not.out1file) &
3171 write (iout,*) "REMD setup"
3173 call card_concat(controlcard,.true.)
3174 call readi(controlcard,"NREP",nrep,3)
3175 call readi(controlcard,"NSTEX",nstex,1000)
3176 call reada(controlcard,"RETMIN",retmin,10.0d0)
3177 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3178 mremdsync=(index(controlcard,'SYNC').gt.0)
3179 call readi(controlcard,"NSYN",i_sync_step,100)
3180 restart1file=(index(controlcard,'REST1FILE').gt.0)
3181 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3182 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3183 if(max_cache_traj_use.gt.max_cache_traj) &
3184 max_cache_traj_use=max_cache_traj
3185 if(me.eq.king.or..not.out1file) then
3186 !d if (traj1file) then
3187 !rc caching is in testing - NTWX is not ignored
3188 !d write (iout,*) "NTWX value is ignored"
3189 !d write (iout,*) " trajectory is stored to one file by master"
3190 !d write (iout,*) " before exchange at NSTEX intervals"
3192 write (iout,*) "NREP= ",nrep
3193 write (iout,*) "NSTEX= ",nstex
3194 write (iout,*) "SYNC= ",mremdsync
3195 write (iout,*) "NSYN= ",i_sync_step
3196 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3199 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3200 if (index(controlcard,'TLIST').gt.0) then
3202 call card_concat(controlcard1,.true.)
3203 read(controlcard1,*) (remd_t(i),i=1,nrep)
3204 if(me.eq.king.or..not.out1file) &
3205 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3208 if (index(controlcard,'MLIST').gt.0) then
3210 call card_concat(controlcard1,.true.)
3211 read(controlcard1,*) (remd_m(i),i=1,nrep)
3212 if(me.eq.king.or..not.out1file) then
3213 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3216 iremd_m_total=iremd_m_total+remd_m(i)
3218 write (iout,*) 'Total number of replicas ',iremd_m_total
3221 if(me.eq.king.or..not.out1file) &
3222 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
3224 end subroutine read_REMDpar
3225 !-----------------------------------------------------------------------------
3226 subroutine read_MDpar
3230 use control_data, only: r_cut,rlamb,out1file
3232 use geometry_data, only: pi
3234 ! implicit real*8 (a-h,o-z)
3235 ! include 'DIMENSIONS'
3236 ! include 'COMMON.IOUNITS'
3237 ! include 'COMMON.TIME1'
3238 ! include 'COMMON.MD'
3241 !el include 'COMMON.LANGEVIN'
3243 !el include 'COMMON.LANGEVIN.lang0'
3245 ! include 'COMMON.INTERACT'
3246 ! include 'COMMON.NAMES'
3247 ! include 'COMMON.GEO'
3248 ! include 'COMMON.SETUP'
3249 ! include 'COMMON.CONTROL'
3250 ! include 'COMMON.SPLITELE'
3251 ! character(len=80) :: ucase
3252 character(len=320) :: controlcard
3257 call card_concat(controlcard,.true.)
3258 call readi(controlcard,"NSTEP",n_timestep,1000000)
3259 call readi(controlcard,"NTWE",ntwe,100)
3260 call readi(controlcard,"NTWX",ntwx,1000)
3261 call reada(controlcard,"DT",d_time,1.0d-1)
3262 call reada(controlcard,"DVMAX",dvmax,2.0d1)
3263 call reada(controlcard,"DAMAX",damax,1.0d1)
3264 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
3265 call readi(controlcard,"LANG",lang,0)
3266 RESPA = index(controlcard,"RESPA") .gt. 0
3267 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
3268 ntime_split0=ntime_split
3269 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
3270 ntime_split0=ntime_split
3271 call reada(controlcard,"R_CUT",r_cut,2.0d0)
3272 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
3273 rest = index(controlcard,"REST").gt.0
3274 tbf = index(controlcard,"TBF").gt.0
3275 usampl = index(controlcard,"USAMPL").gt.0
3276 mdpdb = index(controlcard,"MDPDB").gt.0
3277 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3278 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
3279 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
3280 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
3281 if (count_reset_moment.eq.0) count_reset_moment=1000000000
3282 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
3283 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
3284 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
3285 if (count_reset_vel.eq.0) count_reset_vel=1000000000
3286 large = index(controlcard,"LARGE").gt.0
3287 print_compon = index(controlcard,"PRINT_COMPON").gt.0
3288 rattle = index(controlcard,"RATTLE").gt.0
3289 ! if performing umbrella sampling, fragments constrained are read from the fragment file
3295 if(me.eq.king.or..not.out1file) then
3297 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
3299 write (iout,'(a)') "The units are:"
3300 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
3301 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
3302 " acceleration: angstrom/(48.9 fs)**2"
3303 write (iout,'(a)') "energy: kcal/mol, temperature: K"
3305 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
3306 write (iout,'(a60,f10.5,a)') &
3307 "Initial time step of numerical integration:",d_time,&
3309 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
3311 write (iout,'(2a,i4,a)') &
3312 "A-MTS algorithm used; initial time step for fast-varying",&
3313 " short-range forces split into",ntime_split," steps."
3314 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
3315 r_cut," lambda",rlamb
3317 write (iout,'(2a,f10.5)') &
3318 "Maximum acceleration threshold to reduce the time step",&
3319 "/increase split number:",damax
3320 write (iout,'(2a,f10.5)') &
3321 "Maximum predicted energy drift to reduce the timestep",&
3322 "/increase split number:",edriftmax
3323 write (iout,'(a60,f10.5)') &
3324 "Maximum velocity threshold to reduce velocities:",dvmax
3325 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
3326 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
3327 if (rattle) write (iout,'(a60)') &
3328 "Rattle algorithm used to constrain the virtual bonds"
3332 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
3333 call reada(controlcard,"RWAT",rwat,1.4d0)
3334 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
3335 surfarea=index(controlcard,"SURFAREA").gt.0
3336 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
3337 if(me.eq.king.or..not.out1file)then
3338 write (iout,'(/a,$)') "Langevin dynamics calculation"
3340 write (iout,'(a/)') &
3341 " with direct integration of Langevin equations"
3342 else if (lang.eq.2) then
3343 write (iout,'(a/)') " with TINKER stochasic MD integrator"
3344 else if (lang.eq.3) then
3345 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
3346 else if (lang.eq.4) then
3347 write (iout,'(a/)') " in overdamped mode"
3349 write (iout,'(//a,i5)') &
3350 "=========== ERROR: Unknown Langevin dynamics mode:",lang
3353 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3354 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
3355 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
3356 write (iout,'(a60,f10.5)') &
3357 "Scaling factor of the friction forces:",scal_fric
3358 if (surfarea) write (iout,'(2a,i10,a)') &
3359 "Friction coefficients will be scaled by solvent-accessible",&
3360 " surface area every",reset_fricmat," steps."
3362 ! Calculate friction coefficients and bounds of stochastic forces
3363 eta=6*pi*cPoise*etawat
3364 if(me.eq.king.or..not.out1file) &
3365 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
3367 gamp=scal_fric*(pstok+rwat)*eta
3368 stdfp=dsqrt(2*Rb*t_bath/d_time)
3369 allocate(gamsc(ntyp1),stdfsc(ntyp1)) !(ntyp1)
3371 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
3372 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
3374 if(me.eq.king.or..not.out1file)then
3375 write (iout,'(/2a/)') &
3376 "Radii of site types and friction coefficients and std's of",&
3377 " stochastic forces of fully exposed sites"
3378 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
3380 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),&
3381 gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
3385 if(me.eq.king.or..not.out1file)then
3386 write (iout,'(a)') "Berendsen bath calculation"
3387 write (iout,'(a60,f10.5)') "Temperature:",t_bath
3388 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
3390 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
3391 count_reset_moment," steps"
3393 write (iout,'(a,i10,a)') &
3394 "Velocities will be reset at random every",count_reset_vel,&
3398 if(me.eq.king.or..not.out1file) &
3399 write (iout,'(a31)') "Microcanonical mode calculation"
3401 if(me.eq.king.or..not.out1file)then
3402 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
3404 write(iout,*) "MD running with constraints."
3405 write(iout,*) "Equilibration time ", eq_time, " mtus."
3406 write(iout,*) "Constraining ", nfrag," fragments."
3407 write(iout,*) "Length of each fragment, weight and q0:"
3409 write (iout,*) "Set of restraints #",iset
3411 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
3412 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
3414 write(iout,*) "constraints between ", npair, "fragments."
3415 write(iout,*) "constraint pairs, weights and q0:"
3417 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
3418 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
3420 write(iout,*) "angle constraints within ", nfrag_back,&
3421 "backbone fragments."
3422 write(iout,*) "fragment, weights:"
3424 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
3425 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
3426 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
3429 iset=mod(kolor,nset)+1
3432 if(me.eq.king.or..not.out1file) &
3433 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
3435 end subroutine read_MDpar
3436 !-----------------------------------------------------------------------------
3440 ! implicit real*8 (a-h,o-z)
3441 ! include 'DIMENSIONS'
3442 ! include 'COMMON.MAP'
3443 ! include 'COMMON.IOUNITS'
3444 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
3445 character(len=80) :: mapcard !,ucase
3448 ! real(kind=8) :: var,ene
3451 read (inp,'(a)') mapcard
3452 mapcard=ucase(mapcard)
3453 if (index(mapcard,'PHI').gt.0) then
3455 else if (index(mapcard,'THE').gt.0) then
3457 else if (index(mapcard,'ALP').gt.0) then
3459 else if (index(mapcard,'OME').gt.0) then
3462 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
3463 stop 'Error - illegal variable spec in MAP card.'
3465 call readi (mapcard,'RES1',res1(imap),0)
3466 call readi (mapcard,'RES2',res2(imap),0)
3467 if (res1(imap).eq.0) then
3468 res1(imap)=res2(imap)
3469 else if (res2(imap).eq.0) then
3470 res2(imap)=res1(imap)
3472 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
3473 write (iout,'(a)') &
3474 'Error - illegal definition of variable group in MAP.'
3475 stop 'Error - illegal definition of variable group in MAP.'
3477 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
3478 call reada(mapcard,'TO',ang_to(imap),0.0D0)
3479 call readi(mapcard,'NSTEP',nstep(imap),0)
3480 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
3481 write (iout,'(a)') &
3482 'Illegal boundary and/or step size specification in MAP.'
3483 stop 'Illegal boundary and/or step size specification in MAP.'
3487 end subroutine map_read
3488 !-----------------------------------------------------------------------------
3491 use control_data, only: vdisulf
3493 ! implicit real*8 (a-h,o-z)
3494 ! include 'DIMENSIONS'
3495 ! include 'COMMON.IOUNITS'
3496 ! include 'COMMON.GEO'
3497 ! include 'COMMON.CSA'
3498 ! include 'COMMON.BANK'
3499 ! include 'COMMON.CONTROL'
3500 ! character(len=80) :: ucase
3501 character(len=620) :: mcmcard
3503 ! integer :: ntf,ik,iw_pdb
3504 ! real(kind=8) :: var,ene
3506 call card_concat(mcmcard,.true.)
3508 call readi(mcmcard,'NCONF',nconf,50)
3509 call readi(mcmcard,'NADD',nadd,0)
3510 call readi(mcmcard,'JSTART',jstart,1)
3511 call readi(mcmcard,'JEND',jend,1)
3512 call readi(mcmcard,'NSTMAX',nstmax,500000)
3513 call readi(mcmcard,'N0',n0,1)
3514 call readi(mcmcard,'N1',n1,6)
3515 call readi(mcmcard,'N2',n2,4)
3516 call readi(mcmcard,'N3',n3,0)
3517 call readi(mcmcard,'N4',n4,0)
3518 call readi(mcmcard,'N5',n5,0)
3519 call readi(mcmcard,'N6',n6,10)
3520 call readi(mcmcard,'N7',n7,0)
3521 call readi(mcmcard,'N8',n8,0)
3522 call readi(mcmcard,'N9',n9,0)
3523 call readi(mcmcard,'N14',n14,0)
3524 call readi(mcmcard,'N15',n15,0)
3525 call readi(mcmcard,'N16',n16,0)
3526 call readi(mcmcard,'N17',n17,0)
3527 call readi(mcmcard,'N18',n18,0)
3529 vdisulf=(index(mcmcard,'DYNSS').gt.0)
3531 call readi(mcmcard,'NDIFF',ndiff,2)
3532 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
3533 call readi(mcmcard,'IS1',is1,1)
3534 call readi(mcmcard,'IS2',is2,8)
3535 call readi(mcmcard,'NRAN0',nran0,4)
3536 call readi(mcmcard,'NRAN1',nran1,2)
3537 call readi(mcmcard,'IRR',irr,1)
3538 call readi(mcmcard,'NSEED',nseed,20)
3539 call readi(mcmcard,'NTOTAL',ntotal,10000)
3540 call reada(mcmcard,'CUT1',cut1,2.0d0)
3541 call reada(mcmcard,'CUT2',cut2,5.0d0)
3542 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
3543 call readi(mcmcard,'ICMAX',icmax,3)
3544 call readi(mcmcard,'IRESTART',irestart,0)
3545 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
3548 call reada(mcmcard,'DELE',dele,20.0d0)
3549 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
3550 call readi(mcmcard,'IREF',iref,0)
3551 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
3552 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
3553 call readi(mcmcard,'NCONF_IN',nconf_in,0)
3554 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
3555 write (iout,*) "NCONF_IN",nconf_in
3557 end subroutine csaread
3558 !-----------------------------------------------------------------------------
3562 use control_data, only: MaxMoveType
3565 ! implicit real*8 (a-h,o-z)
3566 ! include 'DIMENSIONS'
3567 ! include 'COMMON.MCM'
3568 ! include 'COMMON.MCE'
3569 ! include 'COMMON.IOUNITS'
3570 ! character(len=80) :: ucase
3571 character(len=320) :: mcmcard
3574 ! real(kind=8) :: var,ene
3576 call card_concat(mcmcard,.true.)
3577 call readi(mcmcard,'MAXACC',maxacc,100)
3578 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
3579 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
3580 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
3581 call readi(mcmcard,'MAXREPM',maxrepm,200)
3582 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
3583 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
3584 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
3585 call reada(mcmcard,'E_UP',e_up,5.0D0)
3586 call reada(mcmcard,'DELTE',delte,0.1D0)
3587 call readi(mcmcard,'NSWEEP',nsweep,5)
3588 call readi(mcmcard,'NSTEPH',nsteph,0)
3589 call readi(mcmcard,'NSTEPC',nstepc,0)
3590 call reada(mcmcard,'TMIN',tmin,298.0D0)
3591 call reada(mcmcard,'TMAX',tmax,298.0D0)
3592 call readi(mcmcard,'NWINDOW',nwindow,0)
3593 call readi(mcmcard,'PRINT_MC',print_mc,0)
3594 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
3595 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
3596 ent_read=(index(mcmcard,'ENT_READ').gt.0)
3597 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
3598 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
3599 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
3600 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
3601 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
3602 if (nwindow.gt.0) then
3603 allocate(winstart(nwindow)) !!el (maxres)
3604 allocate(winend(nwindow)) !!el
3605 allocate(winlen(nwindow)) !!el
3606 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
3608 winlen(i)=winend(i)-winstart(i)+1
3611 if (tmax.lt.tmin) tmax=tmin
3612 if (tmax.eq.tmin) then
3616 if (nstepc.gt.0 .and. nsteph.gt.0) then
3617 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
3618 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
3620 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
3621 ! Probabilities of different move types
3622 sumpro_type(0)=0.0D0
3623 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
3624 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
3625 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
3626 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
3627 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
3628 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
3629 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
3631 print *,'i',i,' sumprotype',sumpro_type(i)
3632 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
3633 print *,'i',i,' sumprotype',sumpro_type(i)
3636 end subroutine mcmread
3637 !-----------------------------------------------------------------------------
3638 subroutine read_minim
3641 ! implicit real*8 (a-h,o-z)
3642 ! include 'DIMENSIONS'
3643 ! include 'COMMON.MINIM'
3644 ! include 'COMMON.IOUNITS'
3645 ! character(len=80) :: ucase
3646 character(len=320) :: minimcard
3648 ! integer :: ntf,ik,iw_pdb
3649 ! real(kind=8) :: var,ene
3651 call card_concat(minimcard,.true.)
3652 call readi(minimcard,'MAXMIN',maxmin,2000)
3653 call readi(minimcard,'MAXFUN',maxfun,5000)
3654 call readi(minimcard,'MINMIN',minmin,maxmin)
3655 call readi(minimcard,'MINFUN',minfun,maxmin)
3656 call reada(minimcard,'TOLF',tolf,1.0D-2)
3657 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
3658 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
3659 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
3660 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
3661 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
3662 'Options in energy minimization:'
3663 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
3664 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
3665 'MinMin:',MinMin,' MinFun:',MinFun,&
3666 ' TolF:',TolF,' RTolF:',RTolF
3668 end subroutine read_minim
3669 !-----------------------------------------------------------------------------
3670 subroutine openunits
3672 use MD_data, only: usampl
3675 use control_data, only:out1file
3676 use control, only: getenv_loc
3677 ! implicit real*8 (a-h,o-z)
3678 ! include 'DIMENSIONS'
3681 character(len=16) :: form,nodename
3682 integer :: nodelen,ierror,npos
3684 ! include 'COMMON.SETUP'
3685 ! include 'COMMON.IOUNITS'
3686 ! include 'COMMON.MD'
3687 ! include 'COMMON.CONTROL'
3688 integer :: lenpre,lenpot,lentmp !,ilen
3690 character(len=3) :: out1file_text !,ucase
3691 character(len=3) :: ll
3694 ! integer :: ntf,ik,iw_pdb
3695 ! real(kind=8) :: var,ene
3697 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
3698 call getenv_loc("PREFIX",prefix)
3700 call getenv_loc("POT",pot)
3701 call getenv_loc("DIRTMP",tmpdir)
3702 call getenv_loc("CURDIR",curdir)
3703 call getenv_loc("OUT1FILE",out1file_text)
3704 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
3705 out1file_text=ucase(out1file_text)
3706 if (out1file_text(1:1).eq."Y") then
3709 out1file=fg_rank.gt.0
3714 if (lentmp.gt.0) then
3715 write (*,'(80(1h!))')
3716 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
3717 write (*,'(80(1h!))')
3718 write (*,*)"All output files will be on node /tmp directory."
3720 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
3721 if (me.eq.king) then
3722 write (*,*) "The master node is ",nodename
3723 else if (fg_rank.eq.0) then
3724 write (*,*) "I am the CG slave node ",nodename
3726 write (*,*) "I am the FG slave node ",nodename
3729 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
3730 lenpre = lentmp+lenpre+1
3732 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
3733 ! Get the names and open the input files
3734 #if defined(WINIFL) || defined(WINPGI)
3735 open(1,file=pref_orig(:ilen(pref_orig))// &
3736 '.inp',status='old',readonly,shared)
3737 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3738 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3739 ! Get parameter filenames and open the parameter files.
3740 call getenv_loc('BONDPAR',bondname)
3741 open (ibond,file=bondname,status='old',readonly,shared)
3742 call getenv_loc('THETPAR',thetname)
3743 open (ithep,file=thetname,status='old',readonly,shared)
3744 call getenv_loc('ROTPAR',rotname)
3745 open (irotam,file=rotname,status='old',readonly,shared)
3746 call getenv_loc('TORPAR',torname)
3747 open (itorp,file=torname,status='old',readonly,shared)
3748 call getenv_loc('TORDPAR',tordname)
3749 open (itordp,file=tordname,status='old',readonly,shared)
3750 call getenv_loc('FOURIER',fouriername)
3751 open (ifourier,file=fouriername,status='old',readonly,shared)
3752 call getenv_loc('ELEPAR',elename)
3753 open (ielep,file=elename,status='old',readonly,shared)
3754 call getenv_loc('SIDEPAR',sidename)
3755 open (isidep,file=sidename,status='old',readonly,shared)
3756 #elif (defined CRAY) || (defined AIX)
3757 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3759 ! print *,"Processor",myrank," opened file 1"
3760 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3761 ! print *,"Processor",myrank," opened file 9"
3762 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3763 ! Get parameter filenames and open the parameter files.
3764 call getenv_loc('BONDPAR',bondname)
3765 open (ibond,file=bondname,status='old',action='read')
3766 ! print *,"Processor",myrank," opened file IBOND"
3767 call getenv_loc('THETPAR',thetname)
3768 open (ithep,file=thetname,status='old',action='read')
3769 ! print *,"Processor",myrank," opened file ITHEP"
3770 call getenv_loc('ROTPAR',rotname)
3771 open (irotam,file=rotname,status='old',action='read')
3772 ! print *,"Processor",myrank," opened file IROTAM"
3773 call getenv_loc('TORPAR',torname)
3774 open (itorp,file=torname,status='old',action='read')
3775 ! print *,"Processor",myrank," opened file ITORP"
3776 call getenv_loc('TORDPAR',tordname)
3777 open (itordp,file=tordname,status='old',action='read')
3778 ! print *,"Processor",myrank," opened file ITORDP"
3779 call getenv_loc('SCCORPAR',sccorname)
3780 open (isccor,file=sccorname,status='old',action='read')
3781 ! print *,"Processor",myrank," opened file ISCCOR"
3782 call getenv_loc('FOURIER',fouriername)
3783 open (ifourier,file=fouriername,status='old',action='read')
3784 ! print *,"Processor",myrank," opened file IFOURIER"
3785 call getenv_loc('ELEPAR',elename)
3786 open (ielep,file=elename,status='old',action='read')
3787 ! print *,"Processor",myrank," opened file IELEP"
3788 call getenv_loc('SIDEPAR',sidename)
3789 open (isidep,file=sidename,status='old',action='read')
3790 call getenv_loc('LIPTRANPAR',liptranname)
3791 open (iliptranpar,file=liptranname,status='old',action='read')
3792 call getenv_loc('TUBEPAR',tubename)
3793 open (itube,file=tubename,status='old',action='read')
3795 ! print *,"Processor",myrank," opened file ISIDEP"
3796 ! print *,"Processor",myrank," opened parameter files"
3798 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
3799 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3800 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3801 ! Get parameter filenames and open the parameter files.
3802 call getenv_loc('BONDPAR',bondname)
3803 open (ibond,file=bondname,status='old')
3804 call getenv_loc('THETPAR',thetname)
3805 open (ithep,file=thetname,status='old')
3806 call getenv_loc('ROTPAR',rotname)
3807 open (irotam,file=rotname,status='old')
3808 call getenv_loc('TORPAR',torname)
3809 open (itorp,file=torname,status='old')
3810 call getenv_loc('TORDPAR',tordname)
3811 open (itordp,file=tordname,status='old')
3812 call getenv_loc('SCCORPAR',sccorname)
3813 open (isccor,file=sccorname,status='old')
3814 call getenv_loc('FOURIER',fouriername)
3815 open (ifourier,file=fouriername,status='old')
3816 call getenv_loc('ELEPAR',elename)
3817 open (ielep,file=elename,status='old')
3818 call getenv_loc('SIDEPAR',sidename)
3819 open (isidep,file=sidename,status='old')
3820 call getenv_loc('LIPTRANPAR',liptranname)
3821 open (iliptranpar,file=liptranname,status='old')
3822 call getenv_loc('TUBEPAR',tubename)
3823 open (itube,file=tubename,status='old')
3825 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
3827 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
3828 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
3829 ! Get parameter filenames and open the parameter files.
3830 call getenv_loc('BONDPAR',bondname)
3831 open (ibond,file=bondname,status='old',action='read')
3832 call getenv_loc('THETPAR',thetname)
3833 open (ithep,file=thetname,status='old',action='read')
3834 call getenv_loc('ROTPAR',rotname)
3835 open (irotam,file=rotname,status='old',action='read')
3836 call getenv_loc('TORPAR',torname)
3837 open (itorp,file=torname,status='old',action='read')
3838 call getenv_loc('TORDPAR',tordname)
3839 open (itordp,file=tordname,status='old',action='read')
3840 call getenv_loc('SCCORPAR',sccorname)
3841 open (isccor,file=sccorname,status='old',action='read')
3843 call getenv_loc('THETPARPDB',thetname_pdb)
3844 print *,"thetname_pdb ",thetname_pdb
3845 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
3846 print *,ithep_pdb," opened"
3848 call getenv_loc('FOURIER',fouriername)
3849 open (ifourier,file=fouriername,status='old',readonly)
3850 call getenv_loc('ELEPAR',elename)
3851 open (ielep,file=elename,status='old',readonly)
3852 call getenv_loc('SIDEPAR',sidename)
3853 open (isidep,file=sidename,status='old',readonly)
3854 call getenv_loc('LIPTRANPAR',liptranname)
3855 open (iliptranpar,file=liptranname,status='old',action='read')
3856 call getenv_loc('TUBEPAR',tubename)
3857 open (itube,file=tubename,status='old',action='read')
3860 call getenv_loc('ROTPARPDB',rotname_pdb)
3861 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
3866 ! 8/9/01 In the newest version SCp interaction constants are read from a file
3867 ! Use -DOLDSCP to use hard-coded constants instead.
3869 call getenv_loc('SCPPAR',scpname)
3870 #if defined(WINIFL) || defined(WINPGI)
3871 open (iscpp,file=scpname,status='old',readonly,shared)
3872 #elif (defined CRAY) || (defined AIX)
3873 open (iscpp,file=scpname,status='old',action='read')
3875 open (iscpp,file=scpname,status='old')
3877 open (iscpp,file=scpname,status='old',action='read')
3880 call getenv_loc('PATTERN',patname)
3881 #if defined(WINIFL) || defined(WINPGI)
3882 open (icbase,file=patname,status='old',readonly,shared)
3883 #elif (defined CRAY) || (defined AIX)
3884 open (icbase,file=patname,status='old',action='read')
3886 open (icbase,file=patname,status='old')
3888 open (icbase,file=patname,status='old',action='read')
3891 ! Open output file only for CG processes
3892 ! print *,"Processor",myrank," fg_rank",fg_rank
3893 if (fg_rank.eq.0) then
3895 if (nodes.eq.1) then
3898 npos = dlog10(dfloat(nodes-1))+1
3900 if (npos.lt.3) npos=3
3901 write (liczba,'(i1)') npos
3902 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
3904 write (liczba,form) me
3905 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
3906 liczba(:ilen(liczba))
3907 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3909 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
3911 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
3912 liczba(:ilen(liczba))//'.mol2'
3913 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3914 liczba(:ilen(liczba))//'.stat'
3916 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
3917 //liczba(:ilen(liczba))//'.stat')
3918 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
3921 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
3922 liczba(:ilen(liczba))//'.const'
3927 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
3928 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
3929 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
3930 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
3931 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
3933 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
3935 rest2name=prefix(:ilen(prefix))//'.rst'
3937 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
3940 #if defined(AIX) || defined(PGI)
3941 if (me.eq.king .or. .not. out1file) &
3942 open(iout,file=outname,status='unknown')
3944 if (fg_rank.gt.0) then
3945 write (liczba,'(i3.3)') myrank/nfgtasks
3946 write (ll,'(bz,i3.3)') fg_rank
3947 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3952 open(igeom,file=intname,status='unknown',position='append')
3953 open(ipdb,file=pdbname,status='unknown')
3954 open(imol2,file=mol2name,status='unknown')
3955 open(istat,file=statname,status='unknown',position='append')
3957 !1out open(iout,file=outname,status='unknown')
3960 if (me.eq.king .or. .not.out1file) &
3961 open(iout,file=outname,status='unknown')
3963 if (fg_rank.gt.0) then
3964 write (liczba,'(i3.3)') myrank/nfgtasks
3965 write (ll,'(bz,i3.3)') fg_rank
3966 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
3971 open(igeom,file=intname,status='unknown',access='append')
3972 open(ipdb,file=pdbname,status='unknown')
3973 open(imol2,file=mol2name,status='unknown')
3974 open(istat,file=statname,status='unknown',access='append')
3976 !1out open(iout,file=outname,status='unknown')
3979 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
3980 csa_seed=prefix(:lenpre)//'.CSA.seed'
3981 csa_history=prefix(:lenpre)//'.CSA.history'
3982 csa_bank=prefix(:lenpre)//'.CSA.bank'
3983 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
3984 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
3985 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
3986 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
3987 csa_int=prefix(:lenpre)//'.int'
3988 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
3989 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
3990 csa_in=prefix(:lenpre)//'.CSA.in'
3991 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
3994 write (iout,'(80(1h-))')
3995 write (iout,'(30x,a)') "FILE ASSIGNMENT"
3996 write (iout,'(80(1h-))')
3997 write (iout,*) "Input file : ",&
3998 pref_orig(:ilen(pref_orig))//'.inp'
3999 write (iout,*) "Output file : ",&
4000 outname(:ilen(outname))
4002 write (iout,*) "Sidechain potential file : ",&
4003 sidename(:ilen(sidename))
4005 write (iout,*) "SCp potential file : ",&
4006 scpname(:ilen(scpname))
4008 write (iout,*) "Electrostatic potential file : ",&
4009 elename(:ilen(elename))
4010 write (iout,*) "Cumulant coefficient file : ",&
4011 fouriername(:ilen(fouriername))
4012 write (iout,*) "Torsional parameter file : ",&
4013 torname(:ilen(torname))
4014 write (iout,*) "Double torsional parameter file : ",&
4015 tordname(:ilen(tordname))
4016 write (iout,*) "SCCOR parameter file : ",&
4017 sccorname(:ilen(sccorname))
4018 write (iout,*) "Bond & inertia constant file : ",&
4019 bondname(:ilen(bondname))
4020 write (iout,*) "Bending parameter file : ",&
4021 thetname(:ilen(thetname))
4022 write (iout,*) "Rotamer parameter file : ",&
4023 rotname(:ilen(rotname))
4026 write (iout,*) "Thetpdb parameter file : ",&
4027 thetname_pdb(:ilen(thetname_pdb))
4030 write (iout,*) "Threading database : ",&
4031 patname(:ilen(patname))
4033 write (iout,*)" DIRTMP : ",&
4035 write (iout,'(80(1h-))')
4038 end subroutine openunits
4039 !-----------------------------------------------------------------------------
4042 use geometry_data, only: nres,dc
4044 ! implicit real*8 (a-h,o-z)
4045 ! include 'DIMENSIONS'
4046 ! include 'COMMON.CHAIN'
4047 ! include 'COMMON.IOUNITS'
4048 ! include 'COMMON.MD'
4051 ! real(kind=8) :: var,ene
4053 open(irest2,file=rest2name,status='unknown')
4054 read(irest2,*) totT,EK,potE,totE,t_bath
4056 ! AL 4/17/17: Now reading d_t(0,:) too
4058 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4061 ! AL 4/17/17: Now reading d_c(0,:) too
4063 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4066 read (irest2,*) iset
4070 end subroutine readrst
4071 !-----------------------------------------------------------------------------
4072 subroutine read_fragments
4076 use control_data, only:out1file
4079 ! implicit real*8 (a-h,o-z)
4080 ! include 'DIMENSIONS'
4084 ! include 'COMMON.SETUP'
4085 ! include 'COMMON.CHAIN'
4086 ! include 'COMMON.IOUNITS'
4087 ! include 'COMMON.MD'
4088 ! include 'COMMON.CONTROL'
4091 ! real(kind=8) :: var,ene
4093 read(inp,*) nset,nfrag,npair,nfrag_back
4095 !el from module energy
4096 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
4097 if(.not.allocated(wfrag_back)) then
4098 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4099 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4101 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4102 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
4104 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4105 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
4108 if(me.eq.king.or..not.out1file) &
4109 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4110 " nfrag_back",nfrag_back
4112 read(inp,*) mset(iset)
4114 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4116 if(me.eq.king.or..not.out1file) &
4117 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4118 ifrag(2,i,iset), qinfrag(i,iset)
4121 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4123 if(me.eq.king.or..not.out1file) &
4124 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
4125 ipair(2,i,iset), qinpair(i,iset)
4128 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4129 wfrag_back(3,i,iset),&
4130 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4131 if(me.eq.king.or..not.out1file) &
4132 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
4133 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
4137 end subroutine read_fragments
4138 !-----------------------------------------------------------------------------
4140 !-----------------------------------------------------------------------------
4144 ! implicit real*8 (a-h,o-z)
4145 ! include 'DIMENSIONS'
4146 ! include 'COMMON.CSA'
4147 ! include 'COMMON.BANK'
4148 ! include 'COMMON.IOUNITS'
4150 ! integer :: ntf,ik,iw_pdb
4151 ! real(kind=8) :: var,ene
4153 open(icsa_in,file=csa_in,status="old",err=100)
4154 read(icsa_in,*) nconf
4155 read(icsa_in,*) jstart,jend
4156 read(icsa_in,*) nstmax
4157 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4158 read(icsa_in,*) nran0,nran1,irr
4159 read(icsa_in,*) nseed
4160 read(icsa_in,*) ntotal,cut1,cut2
4161 read(icsa_in,*) estop
4162 read(icsa_in,*) icmax,irestart
4163 read(icsa_in,*) ntbankm,dele,difcut
4164 read(icsa_in,*) iref,rmscut,pnccut
4165 read(icsa_in,*) ndiff
4172 end subroutine csa_read
4173 !-----------------------------------------------------------------------------
4174 subroutine initial_write
4177 ! implicit real*8 (a-h,o-z)
4178 ! include 'DIMENSIONS'
4179 ! include 'COMMON.CSA'
4180 ! include 'COMMON.BANK'
4181 ! include 'COMMON.IOUNITS'
4183 ! integer :: ntf,ik,iw_pdb
4184 ! real(kind=8) :: var,ene
4186 open(icsa_seed,file=csa_seed,status="unknown")
4187 write(icsa_seed,*) "seed"
4189 #if defined(AIX) || defined(PGI)
4190 open(icsa_history,file=csa_history,status="unknown",&
4193 open(icsa_history,file=csa_history,status="unknown",&
4196 write(icsa_history,*) nconf
4197 write(icsa_history,*) jstart,jend
4198 write(icsa_history,*) nstmax
4199 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4200 write(icsa_history,*) nran0,nran1,irr
4201 write(icsa_history,*) nseed
4202 write(icsa_history,*) ntotal,cut1,cut2
4203 write(icsa_history,*) estop
4204 write(icsa_history,*) icmax,irestart
4205 write(icsa_history,*) ntbankm,dele,difcut
4206 write(icsa_history,*) iref,rmscut,pnccut
4207 write(icsa_history,*) ndiff
4209 write(icsa_history,*)
4212 open(icsa_bank1,file=csa_bank1,status="unknown")
4213 write(icsa_bank1,*) 0
4217 end subroutine initial_write
4218 !-----------------------------------------------------------------------------
4219 subroutine restart_write
4222 ! implicit real*8 (a-h,o-z)
4223 ! include 'DIMENSIONS'
4224 ! include 'COMMON.IOUNITS'
4225 ! include 'COMMON.CSA'
4226 ! include 'COMMON.BANK'
4228 ! integer :: ntf,ik,iw_pdb
4229 ! real(kind=8) :: var,ene
4231 #if defined(AIX) || defined(PGI)
4232 open(icsa_history,file=csa_history,position="append")
4234 open(icsa_history,file=csa_history,access="append")
4236 write(icsa_history,*)
4237 write(icsa_history,*) "This is restart"
4238 write(icsa_history,*)
4239 write(icsa_history,*) nconf
4240 write(icsa_history,*) jstart,jend
4241 write(icsa_history,*) nstmax
4242 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
4243 write(icsa_history,*) nran0,nran1,irr
4244 write(icsa_history,*) nseed
4245 write(icsa_history,*) ntotal,cut1,cut2
4246 write(icsa_history,*) estop
4247 write(icsa_history,*) icmax,irestart
4248 write(icsa_history,*) ntbankm,dele,difcut
4249 write(icsa_history,*) iref,rmscut,pnccut
4250 write(icsa_history,*) ndiff
4251 write(icsa_history,*)
4252 write(icsa_history,*) "irestart is: ", irestart
4254 write(icsa_history,*)
4258 end subroutine restart_write
4259 !-----------------------------------------------------------------------------
4261 !-----------------------------------------------------------------------------
4262 subroutine write_pdb(npdb,titelloc,ee)
4264 ! implicit real*8 (a-h,o-z)
4265 ! include 'DIMENSIONS'
4266 ! include 'COMMON.IOUNITS'
4267 character(len=50) :: titelloc1
4268 character*(*) :: titelloc
4269 character(len=3) :: zahl
4270 character(len=5) :: liczba5
4272 integer :: npdb !,ilen
4276 ! real(kind=8) :: var,ene
4280 if (npdb.lt.1000) then
4281 call numstr(npdb,zahl)
4282 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
4284 if (npdb.lt.10000) then
4285 write(liczba5,'(i1,i4)') 0,npdb
4287 write(liczba5,'(i5)') npdb
4289 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
4291 call pdbout(ee,titelloc1,ipdb)
4294 end subroutine write_pdb
4295 !-----------------------------------------------------------------------------
4297 !-----------------------------------------------------------------------------
4298 subroutine write_thread_summary
4299 ! Thread the sequence through a database of known structures
4300 use control_data, only: refstr
4302 use energy_data, only: n_ene_comp
4304 ! implicit real*8 (a-h,o-z)
4305 ! include 'DIMENSIONS'
4307 use MPI_data !include 'COMMON.INFO'
4310 ! include 'COMMON.CONTROL'
4311 ! include 'COMMON.CHAIN'
4312 ! include 'COMMON.DBASE'
4313 ! include 'COMMON.INTERACT'
4314 ! include 'COMMON.VAR'
4315 ! include 'COMMON.THREAD'
4316 ! include 'COMMON.FFIELD'
4317 ! include 'COMMON.SBRIDGE'
4318 ! include 'COMMON.HEADER'
4319 ! include 'COMMON.NAMES'
4320 ! include 'COMMON.IOUNITS'
4321 ! include 'COMMON.TIME1'
4323 integer,dimension(maxthread) :: ip
4324 real(kind=8),dimension(0:n_ene) :: energia
4326 integer :: i,j,ii,jj,ipj,ik,kk,ist
4327 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
4329 write (iout,'(30x,a/)') &
4330 ' *********** Summary threading statistics ************'
4331 write (iout,'(a)') 'Initial energies:'
4332 write (iout,'(a4,2x,a12,14a14,3a8)') &
4333 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
4334 'RMSnat','NatCONT','NNCONT','RMS'
4335 ! Energy sort patterns
4340 enet=ener(n_ene-1,ip(i))
4343 if (ener(n_ene-1,ip(j)).lt.enet) then
4345 enet=ener(n_ene-1,ip(j))
4357 ist=nres_base(2,ii)+ipatt(2,i)
4359 energia(i)=ener0(kk,i)
4361 etot=ener0(n_ene_comp+1,i)
4362 rmsnat=ener0(n_ene_comp+2,i)
4363 rms=ener0(n_ene_comp+3,i)
4364 frac=ener0(n_ene_comp+4,i)
4365 frac_nn=ener0(n_ene_comp+5,i)
4368 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4369 i,str_nam(ii),ist+1,&
4370 (energia(print_order(kk)),kk=1,nprint_ene),&
4371 etot,rmsnat,frac,frac_nn,rms
4373 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
4374 i,str_nam(ii),ist+1,&
4375 (energia(print_order(kk)),kk=1,nprint_ene),etot
4378 write (iout,'(//a)') 'Final energies:'
4379 write (iout,'(a4,2x,a12,17a14,3a8)') &
4380 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
4381 'RMSnat','NatCONT','NNCONT','RMS'
4385 ist=nres_base(2,ii)+ipatt(2,i)
4387 energia(kk)=ener(kk,ik)
4389 etot=ener(n_ene_comp+1,i)
4390 rmsnat=ener(n_ene_comp+2,i)
4391 rms=ener(n_ene_comp+3,i)
4392 frac=ener(n_ene_comp+4,i)
4393 frac_nn=ener(n_ene_comp+5,i)
4394 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4395 i,str_nam(ii),ist+1,&
4396 (energia(print_order(kk)),kk=1,nprint_ene),&
4397 etot,rmsnat,frac,frac_nn,rms
4399 write (iout,'(/a/)') 'IEXAM array:'
4400 write (iout,'(i5)') nexcl
4402 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
4404 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
4405 'Max. time for threading step ',max_time_for_thread,&
4406 'Average time for threading step: ',ave_time_for_thread
4408 end subroutine write_thread_summary
4409 !-----------------------------------------------------------------------------
4410 subroutine write_stat_thread(ithread,ipattern,ist)
4412 use energy_data, only: n_ene_comp
4414 ! implicit real*8 (a-h,o-z)
4415 ! include "DIMENSIONS"
4416 ! include "COMMON.CONTROL"
4417 ! include "COMMON.IOUNITS"
4418 ! include "COMMON.THREAD"
4419 ! include "COMMON.FFIELD"
4420 ! include "COMMON.DBASE"
4421 ! include "COMMON.NAMES"
4422 real(kind=8),dimension(0:n_ene) :: energia
4424 integer :: ithread,ipattern,ist,i
4425 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
4427 #if defined(AIX) || defined(PGI)
4428 open(istat,file=statname,position='append')
4430 open(istat,file=statname,access='append')
4433 energia(i)=ener(i,ithread)
4435 etot=ener(n_ene_comp+1,ithread)
4436 rmsnat=ener(n_ene_comp+2,ithread)
4437 rms=ener(n_ene_comp+3,ithread)
4438 frac=ener(n_ene_comp+4,ithread)
4439 frac_nn=ener(n_ene_comp+5,ithread)
4440 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
4441 ithread,str_nam(ipattern),ist+1,&
4442 (energia(print_order(i)),i=1,nprint_ene),&
4443 etot,rmsnat,frac,frac_nn,rms
4446 end subroutine write_stat_thread
4447 !-----------------------------------------------------------------------------
4449 !-----------------------------------------------------------------------------
4450 end module io_config